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A METHOD FOR MEASURING (SLOPES OF) THE MASS PROFILES OF DWARF SPHEROIDAL GALAXIES 
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ABSTRACT 

We introduce a method for measuring the slopes of mass profiles within dwarf spheroidal (dSph) galaxies 
directly from stellar spectroscopic data and without adopting a dark matter halo model. Our method combines 
two recent results: 1 ) spherically symmetric, equilibrium Jeans models imply that the product of halflight radius 
and (squared) stellar velocity dispersion provides an estimate of the mass enclosed within the halflight radius 
of a dSph stellar component, and 2) some dSphs have chemo-dynamically distinct stellar sM^components that 
independently trace the same gravitational potential. We devise a statistical method that uses measurements of 
stellar positions, velocities and spectral indices to distinguish two dSph stellar subcomponents and to estimate 
their individual halflight radii and velocity dispersions. For a dSph with two detected stellar subcomponents, 
we obtain estimates of masses enclosed at two discrete points in the same mass profile, immediately defining a 
slope. Applied to published spectroscopic data, our method distinguishes stellar subcomponents in the Fornax 
and Sculptor dSphs, for which we measure slopes T = AlogM/Alogr = 2.6I+Q37 and V = 2.95^39, respec- 
tively. These values are consistent with 'cores' of constant density within the central few-hundred parsecs of 
each galaxy and rule out 'cuspy' Navarro-Frenk- White (NFW) profiles (c/logM/c/logr < 2 at all radii) with 
significance > 96% and > 99%, respectively. Tests with synthetic data indicate that our method tends system- 
atically to overestimate the mass of the inner stellar subcomponent to a greater degree than that of the outer 
stellar subcomponent, and therefore to underestimate the slope T (implying that the stated NFW exclusion 
levels are conservative). 
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1. INTRODUCTION 

Cold dark matter (CDM) halos produced in collisionless 
cosmological A/-body simulations follow a nearly univer- 
sal mass-density profile that diverges toward the center as 
lim r ^n p(r) oc r~ 7 with 7 ^ 1, forming a so-called 'cusp' 
(Dubinski & Carlbergl [19911 INavarro. Frenk & White! [19961 



volving baryons might alter the original structure of a dark 
matter halo (e gjBlumenfhal et alj|1986fc INavarro et aTl l 1996 ; 
El-Zant et all 1200 It iGnedin et al.1 12004t iTonini et all 120061 
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('NF W' hereaf t er), 

iDiemand et all 120051 ISpringel et alj 12008). 



Klypin et al 
Many 

observations aim to test this scenario by using the mea- 
sured motions of dynamical tracers in individual galax- 
ies to constrain slopes of the unde r lying dark mat - 
ter d e nsity profiles (e.g. , iMoord 1 19941 : iFlores & Primac 



1994; de Blok & McGaugh 1997: ISalucci & Burked 1200' 



McGaugh et alJl2001b iSimon et all 120051 Ide Bloldl2010T~ and 
references therein). However, comparisons to cosmological 
models tend to be inconclusive for the simple reason that 
while most cosmological N-body simulations consider only 
dark matter particles, one observes only baryons. Baryons 
complicate not only the measurement of a dark matter den- 
sity profile but also its interpretation within the context of the 
CDM paradigm. Complicating the measurement is the fact 
that any uncertainty (e.g., uncertain stellar mass-to-light ra- 
tios) in the baryonic mass profile propagates to the inferred 
dark matter profile, as the latter is merely the difference be- 
tween dynamical and baryonic mass profiles. Complicating 
the interpretation of even a clean measurement is the possi- 
bility that various poorly-understood dynamical processes in- 
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Governato et alJl2010tlPontzen & Governatoll201 11) . 

One mitigates both of these complications at once by con- 
sidering the Milky Way's (MW's) dwarf spheroidal (dSph) 
satellites, which have the smallest sizes (~ 10 2 " 3 pc), small- 
est baryonic masses (Ly ~ 1Q 3 ~ 7 Lq) and the largest dynami- 
cal mass-to-ligh t ratios (\M /L V ]/\M /Ly]^ > 10) of any ob- 
served galaxies dAaronsonlll983l lMateo 1998, and references 
therein, Gilmore et al. 2007). Dark matter dominates even at 
the centers of dSph gravitational potential wells, implying that 
uncertainties in baryonic mass profiles have negligible im- 
pact on inferred dark matter profiles. Furthermore, dynami- 
cal processes that invoke baryon physics to alter the central 
structure of dark matter halos are subject to strong constraints 
in dSphs, where one finds the smallest baryon densities and 
infers the largest dark m atter densities of any galaxy type 
dPrvor & Kormendvll 19901) . 

Recent work identifies several mechanisms that in prin- 
ciple are capable of altering the central structure — i.e., 
of transforming primal 'cusps' into 'cores' of constant 
density — of cold dark matter halos on dSph-like scales. 
Such mechanisms typically invoke either the dynamical 
coupli ng of the dark matter to rapid baryonic outflows 
(e.g.. iRead & Gilmord 120051 iMashchenko et all 120061 12001 
Ide Souza et alj|201 11) or the transfer of energy and/or angu- 
lar momen tum to the dark matter from massive infalling ob- 
jects (e.g.. ISanchez-Salcedo et al.1 120061 iGoerdt et alJ I2006L 
120 1 Ot lCole et all 120 1 11) . While both channels for baryon- 
driven core creation are physically plausible, their applica- 
tion to real dSphs must eventually satisfy a large body of 
observational constraints. It is therefore instructive to con- 
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sider recent h igh-re solut ion hydrodynamica l simulations by 
ISawala et all (120101) and iParrv et all ( 12011b . These cosmo- 
logical simulations broadly reproduce the observed distribu- 
tions of luminosities, metallicities and stellar kinematics ex- 
hibited by Local Group dSphs. Both groups conclude that the 
baryon-physical processes driving the formation and evolu- 
tion of their simulated dSphs leave the cuspy central structure 
of dSph CDM halos intact. Taking these results at face value, 
it seems then that the Local Group dSphs represent the most 
pristine dark matter halos to which we have observational ac- 
cess. Measurements of the slopes (i.e. 'cusp' versus 'core') 
of dSph mass profiles can therefore provide a uniquely direct 
test of structure formation within the CDM paradigm. 

Pressure-supported stellar components provide the only 
available kinematic tracers in dSphs, but thus far stellar kine- 
matic data have figured only indirectly in co re/cusp inves- 
tigations. For example, iKleyna et al.1 (120031) detect kine- 
matically cold stellar substructure in the Ursa Minor dSph 
and argue that its survival against tidal disruption is more 
likely in a cored as oppose d to a cusped host po tential. 
ISanchez-Salcedo et all d2006l) and lGoerdt eTail J2006) argue 
that the wide spatial distribution of the five globular clusters 
in the Fornax dSph again favors a cored host potential, as 
dynamical friction within a centrally cusped potential would 
have caused the clusters to sink to Fornax's center in less than 
a Hubble time (unless those c lusters had much wider orbits 
initially). On the other hand, Pefiarrubiaet al. (2010) argue 
that the mass-size relation traced by the Milky Way's dSph 
population favors evolutionary scenarios that invoke cusped 
as opposed to cored halofl 

In contrast to the studies mentioned above, here we devise 
a method for measuring the slopes of dSph mass profiles di- 
rectly from stellar spectroscopic data. We proceed by com- 
bining two recent results. First, for a spherically symmetric 
dSph in dynamic equilibrium, the product of halflight radius 
and (squared) velocity dispersion provid es an estimate of the 
mass enclosed w ithin the halflight radius dWalker et alj 2009a; 
iWolf etalJl2Ol0h . Second, some dSphs contain at least two 
chemo-dynamically distinct stellar populations (|Tolstov et alj 
2004; Battaglia eTaTl 120061: iBattaglia e t"alj 120 l ib , each pre- 
sumably tracing the same dark matter potential. Here we 
formulate a mathematical model that uses measurements of 
stellar positions, velocities and spectral indices to distinguish 
two dSph stellar subcomponents and to estimate their individ- 
ual halflight radii and velocity dispersions. For a dSph with 
two detected stellar subcomponents, we obtain estimates of 
masses enclosed at two discrete points in the same mass pro- 
file. Two points define a slope. 

1.1. Stellar Kinematics with Two Numbers 

In principle th e Collisionless Boltzmann Equation (CBE, 
Equation 4.6 of iBinnev & Tremainel [2008) relates the 6- 
dimensional phase-space distribution function, f(r,v), of a 
tracer component to the underlying gravitational potential, 
thereby governing the joint distribution of stellar positions 

4 This result is particularly s ensitive to the masses i nferre d for the Milky 
Way's 'ultrafaint' satellites. McConnachie & Cote (2010) have recently 
shown that the small velocity dispersions observed for many of these systems 
can receive significant contributions from binary orbital motions, a conclu- 
sion supported by the recent direct detec tion of resolved binary motions in the 
Bootes I satellite (Koposov et al. 2011). Downward revision of the intrinsic 
velocity dispersions (and hence masses) of several of the smallest ultrafaint 
dSphs could lead to a size/mass relation for Milky Way satellites that favor s 
cored over cusped dark matter halos (see Figure 1 1 of Penarrubia et al. 2010). 



and velocities for a pressure-supported galaxy in dynamic 
equilibrium. In practice the available dSph data provide in- 
formation in only three dimensions — two spatial dimensions 
orthogonal to the line of sight and one velocity dimension 
along the line of sight. Implementation of the CBE with dSph 
data then requires transformations between 6D and 3D (or 
2D with spherical sym metry) phase-space distributions (e.g., 
IWilkinson et al.ll2002l) . often at significant computational ex- 
pense. 

Many dSph kinematic st u dies (e .g., Wilki nson et al.l [20041 : 
i Strigari et al] 120061 12001 iKochet alJ 120071: lLokasI 120091: 
Walker et al.ll2009ai: [Battaglia et al.ll2008al l2011[) rely instead 
on the Jeans equations, obtained by integrating the CBE over 
velocity space. The spherically symmetric Jeans equation 
specifies the mass profile M(r) — including the contribution 
from any dark matter component — in terms of the stellar den- 
sity profile, v{r), and stell ar velocity dispersion profile, v 2 (r) 
(IBinnev & Tremaindl2008l) : 

1 d - 2-, - 2 GM(r) 
v dr r r L 

where v 2 and v 2 e are components of the velocity dispersion in 
radial and tangential directions, respectively. Confinement of 
dSph stellar velocity data to the component along the line of 
sight leaves the velocity anisotropy — usually quantified by the 
ratio p w \{r) = 1 — Vg(r) jv^if) — poorly constrained, ultimately 
precluding model-independent constraints on the mass profile 
in analyses based on Equation [T] For example, the top panel 
of FigureQ]demonstrates that the projected velocity dispersion 
profile observed for the Fornax dSph can be fit equally well by 
Jeans models that assume either cored or NFW-cusped dark 
matter halos, or if the shape of the dark matter halo is unspec- 
ified, by models that assume the velocity distribution is either 
isotropic, radially anisotropic or tangentially anisotropic. 

The bottom panel of Figure Q] demonstrates that despite 
this well-known degeneracy between mass and anisotropy, the 
various successful Jeans models tend to have the same mass 
enclosed within app r oximately the dSph halflig h t radius (e.g., 
Strigari et al] l2007t iPenarrubia et all l2008af IWalker et alj 
2009at IWolf et al.ll2010l: lAmorisco & Evansll201 lh . Further- 
more, gi ven the general flatn ess of dSph velocity dispersion 
profiles (Walk er et al.ll2007bl) . the value of this mass relates 
simply to the product of velocity dispersion and the halflight 
radius of the adopted surface brightness model. Assuming 
that the stellar component follows a Plummer profile (v(f) oc 
[1 +r 2 /r\]~ 5 l 2 ) and has an isotrop ic (v 2 = Vp) velocity distri- 
bution with constant dispersion, Walker et al. (2009a) derive 
from Equation[T]the simple estimator 

MW^, (2) 

where r/, is the projected halflight radius and a 2 , is the square 
of the velocity dispersion. Equation [2] provides estimates of 
M{rh) that agree well with formal constraints from searches 
of many-dimensional p arameter spaces often invoked in Jeans 
models (see Figure 3 of lWalker et alll200 9a). This agreement 
follows from the fact that for flat velocity dispersion profiles 
and the adoption of a Plummer surface brightness profile, the 
two observables r/, and a 2 , already contain all of the empirical 
information that enters Equation[T] 

Thus Equation [2] provides reasonably model-independent 
estimates of masses enclosed within dSph projected halflight 



A Method for Measuring (Slopes of) the Mass Profiles in dSph Galaxies 



3 



14 
12 



^ 10 r.-^-^j - jl 



Fornax 



Core 

NFW Cusp 

Isotropic (|S onl =0) 

Radial (/S oni = + 0.3) 

Tangential (/S oni =-0.5) 




100 1000 
Projected Radius [pc] 




100 



1000 



[PC] 



FIG. 1. — Top: Projected stellar velocity dispersion profile for the Fornax dSph 
adopted from Walker et al. 1 2009a). Overlaid are spherical Jeans models that assume 
either a cored dark matter halo, an NFW dark matter halo, or if one lets the shape of 
the dark matter halo vary freely, velocity distributions that are either isotropic, radially 
anisotropic, or tangentially anisotropic. Bottom: Enclosed-mass profiles corresponding 
to th e same models. The vertical dotted line indicates Fornax's projected halflight ra- 
dius ( Irwin & HatzidimitrioLi 1995), where the simple estimator specified by Equationl2l 
gives M(r/,) = [5.3 ± 0.9] X 10 Mq , in agreement with the value common to the various 
successful Jeans models. 



radii ( see section l4~.4.1l for an examination of bias). IWolf et aT] 
(2010) show that an alternative estimator, which can be ex- 
pressed as M(|r;,) r; 4rftiTy / G, gives model-independent esti- 
mates of dSph masses at slightly larger radii (approximately 
the deprojected halflight radius) where good-fitting mass pro- 
files intersect with slightly less scatter. For our purposes here, 
all estimators of the form M(^r/,) oc ryiy, where k is some 
constant, are equivalent. 



1.2. Distinct Stellar Subcomponents in dSphs 

The ability to make a robust estimate of the mass enclosed 
within the halflight radius of a dSph stellar population takes 
on additional significance given recent discoveries that at 
least some dSphs have multiple stellar populations, each of 
which can serve as an independent tracer of the underlying 
gravitational potential. Using VLT/FL AMES spectroscopy 
of ~ 400 stars in the Sculptor dSph, iTolstoy et ail (120041) 
find evidence for two ancient stellar subcomponents that fol- 
low distinct distributions in position, velocity and metallic- 
ity. Imposing a metallicity cut to divide their spectroscopic 
sample i nto relatively 'metal-r ich' and 'metal-poor' subcom- 
ponents, Tolsto vet al.l (|2004) find that the metal-rich sub- 
component is more centrally concentrated an d kinematically 
colder than the metal-poo r subcomponent. Battagl ia et al.1 
(2006) and Battagli a et al.l (1201 ll) report similar discoveries 
for the Fornax and Sextans dSphs based on VLT/FLAMES 
samples of ~ 550 and ~ 200 member stars, respectively — in 
both galaxies a relatively metal-rich subcomponent is again 
more centrally concentrated and has smaller velocity disper- 



sion than a meta l-poorer, k i nem atically hotter, more extended 
subcomponent. Ibata et al. (2006) report similar phenomenol- 
ogy based on Keck/Deimos observations of 44 members of 
the CVnl dSp h (although foll ow-up Keck/Deimos observa- 
tions by ISimon & Gehal 12007 do not reproduce this result). 
Several other dSphs, including Carina, Tucana and some M3 1 
satellites display evidence for distinct stellar subcomponents 
in the form of differing s patial distributions of blue- and red- 
horiz ontal branch stars (Harb eck et al.ll200Tt Bellazzi ni et al.l 
120011) . but these separations have not yet been linked to stellar 
kinematics. 

There are two previously-published efforts to model dSphs 
as superpositions of two dynamically independent stel- 
lar subcomponents tracing the same dark matter potential. 
iMcConnachie et aT] (120071) find that the spatial and velocity 
distributions of several dSphs, including satellites of both the 
Milky Way and M31, are compatible with the presence of 
dynamically distin ct stellar subcomponen ts embedded within 
cuspy NFW halos. Batta glia et al.l(l2008al) apply more general 
dark matter halo models, consideri ng cored as well as cuspy 
NFW halos. IBattaglia etail j2008a) find that while both types 
of halos can plausibly host the two stellar subcomponents they 
detect in Sculptor, the cored halo gives a better fit to the falling 
velocity dispersion profile they measure for Sculptor's metal- 
rich subcomponent. 

In contrast to these previous studies, here we do not adopt a 
dark matter halo model. Nor do we use rigid color/metallicity 
cuts t o separate metal-rich and metal-poor subcomponents 
(c.f. IBattaglia et al.l l2008al) . Instead we devise a statisti- 
cal technique for the purpose of separating and quantifying 
the radius, metallicity and velocity distributions followed by 
two distinct dSph stellar subcomponents that independently 
trace the same gravitational potential. Where we can resolve 
a dSph into two distinct stellar subcomponents and measure 
the halflight radii and velocity dispersions of both subcom- 
ponents, we effectively resolve two discrete points in a mass 
profile dominated by dark matter. The slope then follows di- 
rectly. 

2. DATA 

We use the spectro sc opic data published by 
IWalker. Mateo & Olszewski! (120091 W09 hereafter) for 
three of the Milky Way's 'classical' dSph satellites: Carina, 
Fornax and Sculptor. These data were compiled over five 
years (2004-2008) of observations with the Magellan/Clay 
(6.5m) telescope and Michigan/MIKE Fiber Spectrograph 
(MMFS, PI: Mario Mateo, Co-I: Ed Olszewski) at Las 
Campanas Observatory, Chile. MMFS spectra cover the 
range 5140 A - 5180 A at high resolution (R ~ 20000), 
including th e prominent magnesiu m triplet (MgT) absorption 
feature; see IWalker et al.l d2007al W07 hereafter) and W09 
for details of observations and data reduction. The resulting 
data set includes line-of-sight velocities and spectral indices 
measured individually for ~ 6000 red giant candidates in the 
three dSphs. For the present study we exclude the ~ 500 
stars for which W09 measure a velocity but do not report a 
value for the magnesium index (due to insufficient signal). 
Table Q] gives central coor dinates, helio centric distances 
and integrated luminosities (Mateo 1998) for these three 
galaxies and lists the numbers of stars (including foreground 
interlopers) in the MMFS spectroscopic samples used here. 

For three reasons we exclude from the present study the 
available MMFS data for the Sextans dSph. First, the MMFS 
data set for Sextans is relatively small, containing veloci- 
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ties for just ~ 950 stars, of which only ~ 450 are Sextans 
members. Second, the method we introduc e be low requires 
calibrated broadband photometry (Section 12.11 ). and as re- 
ported by W09, poor observing conditions precluded accurate 
calibration of the photometry used to select Sextans targets. 
Third, in order to compensate for this shortcoming, W09 re- 
port calibrated photometry for ~ 430 of their Sextans targets 
(~ 300 members) using the catalog o flLee et all d2003). How- 
ever, this catalog covers only a few percent of the total surface 
area of Sextans, complica ting the correction for sampling bias 
that we discuss in Section l2~3l Thus a comparable analysis of 
Sextans, which we leave for future work, would require more 
complete photometry and/or careful accounting for these de- 
ficiencies. 

2. 1 . Velocities and Mg- Triplet Indices 

The MMFS velocity sample has a median measurement er- 
ror of ~ 2 km s , smaller than the internal stellar veloc- 
ity dispersions (oy > 6 km s" 1 ) measured for these three 
dSphs. Although their observations were designed primar- 
ily to measure velocities, W07 demonstrate that spectral 
indices — i.e., pseudo-equivalent widths of resolved Fe and 
Mg absorption lines — measured from MMFS spectra corre- 
late with stellar atmospheric parameters such as effective tem- 
perature, surface gravity and metallicity (see Figures 19-20 
of W07). W07 combine indices measured separately for 
magnesium-triplet lines at 5167 A and 5173 A into a com- 
posite MgT index, denoted SMg, that is qualitatively simi- 
lar to the composite calcium-triplet index, ECa, often used 
to infer stell ar metallicities from near-infrared calcium-triplet 
spectra (e.g.. lArmandroff & Da Costal 1 99 UlKoch et alj2006t 
B atta glia et alJl2008bH Starkenbur g et alJl2010h . 

Like the SCa index, the EMg index provides a measure 
of stellar-atmospheric metal abundance, provided that empir- 
ical calibrations adequately remove the dependences of Mg 
opacity on effective temperature and surface gravity. For red 
giant stars that brighten and cool as they expand, these depen- 
dences translate (for stars of a given metallicity) into a nearly 
linear empirical relationship between SMg and luminosity. 
W09 (their Figure 3) quantify this relationship using MMFS 
observations of red giants in six globular clusters that span 
the metallicity range -2.0 < [Fe/H] < -0.5 and have negligi- 
ble (for our purposes' ) internal metallici ty spr eads. Following 
iRutledge et afl (fl997h and iKoch et all d2006h . W09 simulta- 
neously fit the data for each cluster with straight lines sharing 
a common slope, obtaining 



SMg = -(0.079 ± 0.002)(V- V HB ) + SMg' , 



(3) 



where V — Vhb is the offset in V-band luminosity from the hor- 
izontal branch. The slope quantifies the dependence of opac- 
ity on effective temperature and surface gravity, using lumi- 
nosity as a proxy. The intercept, or 'reduced' index EMg' — 
henceforth denoted W' — represents the value of SMg that the 
star would have if it had the surface gravity and temperature 
of a horizontal branch star. Then taking the empirical calibra- 
tion given by Equation [3] at face value, red giants of similar 
metallicity should have similar W'. 

Using the horizontal-branch magnitudes listed in Table Q] 
we apply Equation [3] to obtain reduced magnesium indices 
W' for all dSph stars in the Magellan/MMFS sample. Since 
we are concerned with stellar-atmospheric chemistry only as 
a diagnostic with which to distinguish stellar subcomponents 
independently of their velocity distributions, we do not re- 



quire further calibration of W' values onto an absolute metal- 
licity scale. In what follows we shall use W' as an indicator 
of relative metallicity and we shall refer to subcomponents 
distinguished by W' as 'metal-rich' and 'metal-poor'. 

2.2. Membership 

W07 (see their Figure 1) chose dSph stellar targets from se- 
lection boxes overlaid on red giant branches (RGBs) apparent 
from V , / photometry of each system. These selection regions 
include contamination from foreground Milky Way stars (typ- 
ically late-type dwarfs) and background point sources (unre- 
solved galaxies, quasars). Fortunately, bona fide dSph mem- 
bers follow conspicuous distributions in velocity, spectral in- 
dex and position, which helps to distinguish them from con- 
taminant populations. 

The MMFS data published by W09 include for each 
star a probability of dSph membership that is derived from 
an exp ectation-maximization (EM) algorithm (Walk er et al.1 
2009b), similar to algorithms pre viously used to determine 
membership in o pen clusters (e.g..lSandersll971l) . The EM al- 
gorithm used by IWalker et al.l d2009bl) adopts the foreground 
velocity distributio n given by the Besancon Milky Way model 
(Rob in etal.1 12003) and assumes a single dSph stellar com- 
ponent. Because the latter assumption obviously is at odds 
with the two-subcomponent models we shall consider here, 
we use W09's membership probabilities only to provide non- 
parametric estimates of position, velocity and Mg-index dis- 
tributions followed by non-members (section 13.5b . The ac- 
curacy with which we recove r the input member fraction in 
tests of our method (Section 14.41 ) indicates that these esti- 
mates are not unduly influenced by the subtle differences be- 
tween foreground probabilities derived from single and multi- 
component dSph models. 

2.3. Spatial Selection Bias 

Finally, in light of the fact that here we use the positions of 
stars in the MMFS samples to estimate halflight radii of the 
stellar populations to which they belong, we must consider 
the fact that the spatial distributions of MMFS-observed stars 
may differ from those of the populations from which they are 
drawn. Each MMFS configuration is limited to < 256 targets 
within a 25' field of view, and therefore the distribution of 
stellar positions is not sampled randomly. 

In order to compensate for selection effects due to the in- 
evitable peculiarities of spectroscopic spatial sampling, we as- 
sign a selection probability to each star according to the local 
ratio of the number of observed stars, dN b s , to the number of 
target candidates, dN CW i, selected according to W07's photo- 
metric criteria (the latter number includes both observed and 
unobserved stars). We estimate this ratio as a function of pro- 
jected radius by smoothing the data and target catalogs with 
Gaussian kernels: 



JVob., 



^exp 



w(R) = 



i=\ 



dN ohs (R) 



dN CRnd (R) ^ 

^exp 



1 (Rj-R) 2 
"2 k\ 



1 (Ri-R) 2 
~2 kj 



(4) 



where the hat ( " ) symbol denotes a quantity estimated via 
kernel smoothing. Figure [2] displays w(R) curves for each 
dSph and for possible choices of bandwidth over the range 
0.1 < fci/arcmin) < 10 — these smoothing scales are smaller 
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TABLE 1 

Positions, luminosities and numbers of Magellan/MMFS-observed 
stars for southern classical dsphs* 



dSph RA (J2000) Dec. (J2000) Distance My M-.mi: ^sample 
[hh:mm:ss] [dd:mm:ss] [kpc] [mag] [mag] 



Carina 06:41:37 -50:58:00 101 ±5 -9.3 20.9 1481 

Fornax 02:39:59 -34:27:00 138 ±8 -13.2 21.3 2603 

Sculptor 01:00:09 -33:42:30 79±4 -11.1 20.1 1497 



Central coordinates, distances and absolute magnitudes are adopted from the review of Mateo 

Jp9§. * 

references for horizontal-branch magnitudes: Carina IKoch et aU2006l) , For nax and Sculptor 
llrwin & Hatzidimitriou 1995; Battaglia et al. 2008a) 



than the scale radius of the composite stellar component and 
larger than the (projected) mean free path of sampled stars. 
For the present work we adopt k\ = 2 arcmin (dotted red 
curves in Figure |2J, but we have confirmed that our results 
and conclusions are not qualitatively sensitive to this choice. 

3. METHOD 

We model each dSph mathematically as the superposition 
of two chemo-dynamically distinct stellar subcomponents ob- 
served through foreground contamination. For 'metal-rich' 
and 'metal-poor' stellar subcomponents we adopt simple 
parametric models to describe the distributions of projected 
radius, line-of-sight velocity and magnesium index (Sections 
13.21 - 13 -4b . For the foreground contamination we estimate 
these distributions non-parametrically by smoothing the spec- 
troscopic data according to the publ ishe d probability of (non-) 
membership for each star (Section [3.5b . We constrain model 
parameters using a s tandard Markov-Chain Monte Carlo al- 
gorithm (Section [i3~7l i. The algorithm returns estimates of pa- 
rameters including halfiight radii and velocity dispersions for 
the two stellar subcomponents, which simultaneously provide 
estimates of the mass enclosed at the halfiight radius of each 
subcomponent (Equation |2|i as well as the slope of the mass 
profile: 



AlogM log[M(r h , 2 )/M(r M )] 



Alogr 



l0g|>h,2/Vh,l] 



1 + 



log[/h,2Ah,l 



(5) 

The last expression on the right-hand side follows from 
Equation [2] or, more generally, from any mass estimator 
of the form M(^r^) oc r^a 2 where k is some constant (e.g . 
Penarrubia et al.l2008atlWalker et al.l2009atlWolf et al.l2010b 
Amorisco & Evans 201 1), and makes explicit the fact that the 
only physical quantities relevant to our estimate of the slope 
T are the sizes and velocity dispersions of the two stellar sub- 
components. 

3.1. Likelihood Function 

We require a mathematical model that will let us distin- 
guish and quantify the properties of two independent stel- 
lar subcomponents in the same dSph. Suppose /?[(/?, V,W') 
and p2(R,V,W') describe joint probability distribution^ of 

5 These probability densities are statistical distribution functions de- 
fined such that, for example, the fraction of subcomponent- 1 stars that 
have position in the interval R,R + dR, velocity in the interval V,V + dV 
and reduced magnesium index in the interval W' ,W' + dW' is given by 
p\(R,V,W')dRdVdW'. 



projected radius, line-of-sight velocity and reduced magne- 
sium index for metal-rich and metal-poor dSph subcompo- 
nents, respectively. Further suppose that pmw(R, V,W') is the 
joint probability distribution followed by foreground Milky 
Way stars that satisfy the color/magnitude criteria used to se- 
lect dSph targets. Our data set {R h Vi,W-}^' samples all 

three stellar populations. Let the vector S represent a set 
of free parameters that specifies models for p\(R,V,W') and 
P2(R, V, W) as well as the fractions f\ = Ni /(Nj +N2 +N M w) 
and fi = N2/ (N\ +N2 +Nmw) of stars in metal-rich and metal- 
poor st ellar subcomponents, respectively. Recalling from 
Section 12.31 that w(R) indicates the probability that an RGB 
candidate at radius R is actually observed, then given a model 
specified by S the data set has likelihood 



L({^,y i ,^'}^r lc |5): 

Sample 

/l 

+/2 



n 



w{Rdpi{R h Vi,W[) 



+(1-/1-/2) 



/// w(R)p l (R,V,W')dRdVdW 

w(Ri)p 2 (Ri,Vi,W/) 
J J J w(R)p z (R,V,W')dRdVdW 
w{Rdpyw{RuVi,W!) 



JJJ w(R)p M w(R,V,W')dRdVdW' 



(6) 



The normalizing constant in the denominator of each term 
provid es a weight that compensates for the radial samplin g 
bias dGill et alj!98g;IWanp et alT20lBtlMart"inez et alj|201lb . 

In the following subsections we specify mathematical mod- 
els for p\, P2 and Pmw that l et us evaluate the overall likeli- 
hood given by Equation[6] In selecting from what is in princi- 
ple an unlimited number of possible models, we opt for sim- 
plicity over elegance. That is, in order to quantify distribu- 
tions of positions, velocities and spectral indices we choose 
simple mathematical expressions that let us specify the like- 
lihood function analytically and without introducing an un- 
wieldy number of free parameters. Since one can think of 
p\, p2 and pmw as chemo-dynamical distribution functions, it 
is important to realize that the mathematical models that we 
adopt do not necessarily satisfy the Collisionless Boltzmann 
Equation and thus do not necessarily correspond to physical 
models. It will therefore become necessary to explore the sys- 
tematic errors introduced by the particular mathematical mod- 
els that we adopt (see Section|4]i. 

3.2. Position Distributions for dSph stellar subcomponents 
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FIG. 2. — Spatial sampling bias in the Magellan dSph spectroscopic samples (Walker, Mateo & Olszewski 2009) adopted here. Panels plot the observational selection probability, 
w = dN bsJ dNcsnd, as a function of projected radius, estimated via kernel smoothing (EquationRV Different linestyles correspond to different smoothing bandwidths. For this study we 
adopt the w(R) estimates that correspond to bandwidth k\ = 2 arcmin (dotted red); our results and conclusions are not sensitive to this choice. 



dSph stellar surface densit y (number of star s per un i t area) 
profile s typic a lly are fit by [Plummer (I1911I). iKind (ll962i) 



Irwin & Hatzidimitriou 
Belokurov et all 12007b 



and/o r ISersid (11968I). profiles (e.g., 
199 5I; iMcConnachie & Irwiiil l2006t 
Martin et al.l2008l) . For simplicity we model the spatial distri- 
butions of member subcomponents with spherically symmet- 
ric Plummer profiles, for which the only free parameter is the 
projected halflight radius^, r/,. Then, for example, the metal- 
rich subcomponent (hereafter denoted with the subscript 1) 
has stellar surface density given by 



Si(/f)=- 



'I, A 



[l+R 2 /, 



(7) 



MJ 



(and similarly for the metal-poor subcomponent, hereafter 
denoted with the subscript 2), where N\ is the num- 
ber of stars in subcomponent 1. For a spherical system 
with stellar surface density profile T,(R), the 1-D proba- 
bility distribution of projected radii is given by pu(R) = 

d/dR[f*T,(S)SdS/ f™T,(S)SdS], which for the Plummer 
profile becomes 



Pra(R) = 



(l+R 2 /i 



(8) 



h.\) 



(and similarly for the metal-poor subcomponent). 



3.3. Velocity Distributions for dSph stellar subcomponents 

Again for simplicity, we assume that a given dSph stellar 
subcomponent has an intrinsic line-of-sight velocity distribu- 
tion that is Gaussian and independent of radius. Then the 1-D 
probability distribution of velocities is given by 



Pv,i(V,a*,6*) = 



1 



2n(<T 2 j 



+4) 



: exp 



nv- 

~2 



-(V) a , :S ,) 2 



'v.\ 



+4 



(9) 

(and similarly for the metal -poor subcomponent), where ay.i 
is the intrinsic velocity dispersion and ey is the velocity 
measurement error. The dependence of the mean velocity 
(y)a,.s, on position follows from the fact that the velocity 
data are given in the heliocentric rest frame (HRF, i.e., the 

6 The projected halflight radius is the radius of the circle enclosing half of 
the stars as observed in projection on the plane of the sky. 



rest frame of an observer instantaneously comoving with the 
Sun). Conversion to the dSph rest frame depends on relative 
motion between Sun and dSph, the transverse components of 
which are difficult to measure even from Hubble Space Tele- 
scope (HST) images sep arated by several years (Piatek et al. 
Hj0l [2003, 20061 l2Q07h . Furthermore, the line-of-sight ve- 
locity offset between HRF and the dSph rest frame varies sys- 
tematically along the vector of any transverse motion between 
Sun and dSph, as this motion projects more strongly along 
lines of sight that are farther from the dSph center. 

In order to ac c ount for this 'perspective effect' (e.g., 
[ Feast et all [19611: Ivan der Marel et al] l2002h we follow 
iKaplinghat & Strigaril (12008) and Walke r et akfeoOSl) in con- 
sidering a dSph at heliocentric distance D with central equa- 
torial coordinates (ao,8o), HRF line-of-sight velocity Vjj and 
HRF proper motion (fi a ,fis). At the location of a star at 
(a*, <5*) the dSph's syste mic HRF velocity along the line of 
sight is (see Appendix of Walker et al. 2008 for a derivation) 

(V) a „s, = 

cos <5* sin a* (Vo cos So sin ao + D/i a cos So cos ao 
-D/is sin So sin ao) + cos <5* cos a* (Vo cos So cos ao 
—D/is sin So cos ao — D[i a cos So sin ao) 

+ sin S* (Vo sin So + D[i$ cos <$o)(.l 0) 

3.4. Mg-Index Distributions 

We assume that the 1-D distributions of reduced magne- 
sium index W' in metal-rich and metal-poor dSph subcompo- 
nents are Gaussian and independent of radius such that, for 
example, the probability distribution of magnesium indices 
for the metal-rich subcomponent is 



Pwa(W') = 



1 



2lT ( (J W'A +e W 



: exp 



1 (W 

~2 



■OOi) 2 



a 2 +f 2 



(ID 



(and similarly for the metal-poor subcomponent), where 
(W')i is the mean spectral index, ay/'A is the intrinsic spectral 
index dispersion and ew' is the measurement error. 

3.5. Foreground Distributions 

Compared to dSph member populations, the Milky Way 
foreground has a more uniform spatial distribution, a wider 
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velocity distribution, and systematically larger magnesium in- 
dices. We quantify each of these distributions directly from 
the MMFS data using estimates of membership probabilities, 
Pmem, tabulated by W09. We then smooth the data with Gaus- 
sian kernels to estimate probability distributions of radius, ve- 
locity and reduced magnesium index for the foreground con- 
tamination present in the data: 



^sample 

E 



/>«,Mw(K) - 



2nkl 



-exp 



1 (R~Rf 
~2 k\ 



^sample 

^ ^ (1 Rmem,i) 
i=l 



Sample _ 

\y K J- 1 mem./ 



PV.Mw(Y) - ' 



exp 



2ire 



v.i 



1 (Vi-Vf 



1 V.i 



Sample 

^ (1 Rmem,i) 



PWMwiW) = ■ 



A'samnlc 
_ p ' 1 _ p 

^ s 1 1 mem./ 

, = i J2tt€ 2 w , a 



exp 



1 (W/-W) 2 



2 el, 



-, (12) 



Sample 

^ y (1 ^mem,f) 

!=1 

where ey and e^" are again measurement errors. In general the 
distributions of foreground stars do not change significantly 
over the region (~ 1 square degree) subtended by a dSph, so 
the only spatial scales built into the data are those correspond- 
ing to the dSph struc tural parameters and observational selec- 
tion effects (Section [2.31 ). In order to avoid introducing ad- 
ditional scales into our analysis, we set the spatial smoothing 
bandwidth, k 2 , equal to the halflight radius measured for the 
composite d Sph stellar component. We ado pt halflight radii 
measur ed by I rwin & Hat zidimitrioul d 19951) and tabulated in 
the erratum to Walk er et alj (|2009a). For our tests with syn- 
thetic data (Section [4]i we adopt the halflight radius that we 
estimate from the composite stellar population before consid- 
ering two-subcomponent models. 

3.6. Likelihood Function Revisited 

We have now specified parametric models for the ID prob- 
ability distributions of position, velocity and magnesium in- 
dex for both dSph stellar subcomponents and the Milky Way 
foreground contamination. Because our adopted models for 
velocity and magnesium index distributions are independent 
of radius, the joint probability distributions are sepa rable into 
products of the 1-D distributions adopted in Sections [3~2l - l3.51 
such that 

P 1 (r, v, W) = PR^pv.mpwA (w'y, 

P2(R ) V,W') = p R>2 (R)pv,2(y)p w > >z (W'); 

Pmw(R, V,W') = Prmw(R)Pv,mw(V)pw,mw( w ')- ( 13 ) 

After making these substitutions and recognizing 
that / p v .i(V)dV = J p v , 2 (V)dV = J p WA (W')dW' = 
J Pw,2(W)dW' = 1 by construction, the likelihood given by 



Equation|6]becomes 
L({R h V h W;}^\S)-- 

AU,plc r 



n 



fx 



w{Ri)PRj{Ri)PvAVi)Pw,i{W!) 

/ °° w(R)p RA (R)dR 
w(Ri)p R ,2(Rdpv,2(VdPw',2(W!) 



!™ w{R)p R , 2 (R)dR 



+(1 ~fl -f2)PMW.R(Ri)PMW,v(Vi)pMW.W'(Wi') 



(14) 



For the term representing the Milky Way fore- 
ground we have substituted the probability distri- 
butions estimated by smoothing the data (section 
13.51 ), which necessarily include the effects of sam- 
pling bias such that Pmw,j?(^/)Pmw,vW)Pmw,w"(W/) 

w(^)pMwW-,V,,W/)[/ / Jw(R) P Mw(R,V,W')dRdVdW'Y 1 . 

3.7. Markov-Chain Monte Carlo Technique 

We require a total of twelve free parameters in order to eval- 
uate the likelihood given by Equation [14] Eight parameters 
specify for both member subcomponents the halflight radii 
(Equation [8]), velocity dispersions (Equation [9]), and means 
and variances of the reduced Mg-index distributions (Equa- 
tionHHi. Two parameters, / mem = (M +N 2 )/(Ni +N 2 +N M w) 
and / sub = (Ni/(Ni +N 2 ), specify the fractions f x = f mem f sub 
and /2 = /mem(l _ /sub)- The final two parameters, /i Q and p,s, 
specify the two components of the dSph proper motion (Equa- 
tion[lO]). 

For all parameters we adopt uniform priors over ranges that 
include all reasonable values (Table [2]). Notice that we spec- 
ify halflight radii of the stellar subcomponents in terms of 
free parameters log 10 [r/,.2/pc] and the ratio [r^i/r^L and 
we specify the mean reduced Mg indices in terms of free 
parameters (W)i/A and the difference ((W')i - (W') 2 )/k. 
This formulation lets us apply priors < [//,, 1/77, 2] < 1 and 

< ( ( W } 1 - ( W ) 2 ) / A < 3 . In other words, we assume that the 
metal-rich subcomponent is more centrally concentrated than 
the metal-poor s ubcomponent, as in dicated by previous stud- 
ies of Sculptor (Tol stoy et alj 12004)). Fo rnax (Battagl ia et alj 
2006) and Sextans (Battagli a et alJl20Tlb . 

We sample the large parameter space using standard 
Markov-Chain Monte Carlo techniques . Specifically we use 
the Metropolis -Hastings algorithm (Metr opolis et alJ 11953b 
lHastings| [l970). which samples the parameter space accord- 
ing to the following prescription: i) from the current loca- 
tion in parameter space S„, draw a prospective new loca- 
tion, S', from a Gaussian proposal density centered on S„; 
ii) evaluate the ratio of likelihoods at S„ and S'; and iii) if 
L(S')/L(S„) > 1, accept such that S n+ i = S', else accept with 
probability L(S')/L(S n ), such that S n+ \ = S' with probability 
L(S')/L(S„) and S n+l = S„ with probability 1 -L(S')/L(S n ). 

We implement the Metropolis-Hastings algorithm using th e 
adaptive MCMC engine CosmoMCQ dLewis & Bridldl2002l) . 
CosmoMC provides a generic sampler that periodically up- 
dates the proposal density according to parameter covariances 
in order to optimize the acceptance rate. For a given galaxy 
we run four chains simultaneously, stopping when either the 

7 available at http://cosmologist.info/cosmomc 
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TABLE 2 

Boundaries of Uniform Priors for Free Parameters in Likelihood Function 



Parameter 


Minimum 


Maxium 


Description 


/mem 





1 


= (Ni +A?2)/(JVi +N2 +N MW ), fraction of stars belonging to dSph 


/sub 





1 


= N[ 1 (N\ +N2), fraction of members belonging to MR subcomponent 







1 


ratio of halflight radii for metal-rich (MR) and metal-poor (MP) subcomponents 


l°glC)[''''.2/P c ] 


u 


A S 
■4. J 


halflight radius of MP subcomponent 


(W'),/A 


-3 


+3 


mean reduced Mg index of MR subcomponent 







3 


offset of mean Mg indices 


logio[°"w',iA 2 ] 


-5 


+ 1 


squared dispersion of reduced Mg index, MR subcomponent 




-5 


+ 1 


squared dispersion of reduced Mg index, MP subcomponent 


log 10 [4,7(kmV 2 )] 


-5 


+5 


squared velocity dispersion, MR subcomponent 


log 10 [(7^ 2 /(kmV" 2 )] 


-5 


+5 


squared velocity dispersion, MP subcomponent 


fi a / (mas / century) 


-1000 


+1000 


RA proper motion of dSph 


/j,5 /(mas /century) 


-1000 


+1000 


Dec. proper motion of dSph 



variances of parameter values across the four chains become 
less than 1% of the mean of the variances. Our chains typi- 
cally require of order ~ 10 5 steps to satisfy this convergence 
criterion. 

We process the chain output in two ways. First, in or- 
der to allow for a 'burn-in' period during which the chains 
evolve from initial values to regions of high likelihood, we 
discard the first half of accepted points. Second, because the 
nature of the algorithm introduces correlations between adja- 
cent accepted points, we 'thin' the chains by passing only ev- 
ery 25 th accepted point to a single 'final' chain. The MCMC 
method is designed such that points in the final chain ran- 
domly sample the posterior probability distribution function 
(PDF) in the full twelve-dimensional parameter space. One 
obtains marginalized probability distributions for any single 
parameter, or combination of parameters, simply by counting 
the number of points in the final chain that fall wi thin binned 
ranges of parameter values dLewis & B ridle 2003). 

The final chains sample, among other distributions, the one- 
dimensional posterior PDFs of free parameters log 10 [r/,,2/pc], 
[fh,\/r h ,2], log 10 [cTyj/(km 2 s- 2 )] and log 10 [^ 2 /(kmV 2 )]. 
For every point in the final chains we apply Equation [5] to 
obtain the corresponding slope of the mass profile. The distri- 
bution of these slopes then represents our constraint on T. 

4. TESTS 

Our adoption of Plummer profiles and Gaussian distribu- 
tions to characterize stellar positions, velocities and spectral 
indices (Section [3]) is motivated by preferences for simplic- 
ity and mathematical convenience. While they are frequently 
invoked in analyses of dSph data, there is of course no guaran- 
tee that these particular models provide accurate descriptions 
of real dSphs. In fact, jointly these models can correspond to 
chemodynamical distribution functions that are unphysical. 

For example, our method twice employs the assumption 
that both dSph stellar subcomponents have constant velocity 
dispersion. This assumption is implicit i n the mass estimator 
given by Equation [2] (see Section 3.6 of Wa lker et al.i r2009a) 
and explicit in the models adopted for subcomponent velocity 
distributions (Equation [9]). While the composite stellar popu- 
lations of dSphs generally exhibit flat velocity dispersion pro- 
files over the range of radii where data are available, all (New- 
tonian) equilibrium dynamical models require that tracer ve- 
locity dispersion profiles eventually decline at large radius. 



Furthermore, even in the regions where data indicate that the 
composite dSph stellar populations have flat velocity disper- 
sion profiles, there is no guarantee that the velocity dispersion 
profiles of eit her/both stellar subcomponents are flat individ- 
ually. Indeed, Batta glia et alJ (l2008al) find that the metal-rich 
subcomponent of the Sculptor dSph has a falling velocity dis- 
persion profile (see their Figure 3). 

We must therefore examine the reliability of our method 
when it operates on data sets representing realistic dynami- 
cal systems that violate our simplistic modeling assumptions. 
In order to accomplish this task, we apply our method to a 
series of synthetic data sets and compare the resulting esti- 
mates of stellar population parameters, masses and slopes to 
known input values. We construct a given synthetic data set 
by sampling the superposition of three (two dSph-like mem- 
ber subcomponents plus foreground contamination) chemo- 
dynamical distribution functions. The phase-space distribu- 
tion functions of both member subcomponents correspond to 
independent (except for the obvious constraint that both have 
the same gravitational potential) physical dynamical models 
in which the stellar subcomponent traces a gravitational po- 
tential that is generated by a dark matter halo. Here we de- 
scribe these dynamical models, the generation of synthetic 
data sets, and the systematic behaviors of the errors we iden- 
tify by applying our method to these synthetic data. 

4.1. Dynamical Models for Individual Stellar Subcomponents 

In order to generate test cases we consider physical dynami- 
cal models in which two stellar subcomponents independently 
trace the same dark matter potential. We build these models 
simply by combining the two dynamical models that describe 
the stellar subcomponents individually. We consider individ- 
ual stellar subcomponents that are each dist ributed accordin 
to a genera lized Hernquist density profile (Hernquis3[l99~ 
lZhaolll996h . 



v*(r) = !/()[ — 



r 


1+ (fr 




' * 



(7„-/3„)/a. 



(15) 



and we consider dark matter halos with density profiles that 
take the same form, 



PdmW = po ( 

, '"DM 



(7dm ~/3dm)/ »dm 



(16) 
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These profiles have independent parameters specifying nor- 
malization, scale radius, inner logarithmic slope (7, subscripts 
omitted for brevity), outer logarithmic slope (/?), and the 
sharpness (a) of the transition between the two slopes. 

We aim to test specifically whether our method can distin- 
guish dark matter halos having constant-density cores (7dm = 
0) from those that have NFW-like cusps (7dm =1)- Therefore 
we consider dynamical models in which the central slope of 
the dark matter density profile takes values of either 7dm = 
or 7dm = 1. We hold fixed other halo parameters at scale ra- 
dius tdm = 1 kpc, outer slope /?dm = 3 and «dm = 1 ■ For the 
stars, we consider stellar subcomponents that have structural 
parameters a* } \ = a*,2 = 2 and 7*1 = 7, 2 = 0.1, a range of 
outer slopes /?» =4,5,6, and a range of scale radii r*/roM = 
0.1,0.25,0.5, 1, 1.5 corresponding to various degrees of 'em- 
beddedness' within the dark matter halo. 

In order to generate synthetic data sets we must first calcu- 
late phase-space distribution functions for each stellar sub- 
component in each dark matter halo. For this purpose we 
consider the famil y of spher i cal, an isot ropic dis t ributio n func- 
tions discussed by Osipkovl (119791) and lMerritj dl985l) . These 
models have velocity distributions with anisotropy profiles of 

the form /3ani( r ) = 1 — Vg/v$ = r 2 /(r 2 + rl). We consider values 
for the anisotropy radius r a that give the stellar subcompo- 
nent a velocity distribution that either is isotropic at all radii 
(r„ = 00) or gradually changes from isotropic at small radii 
to radially anisotropic at large radii (r„ = r»). Having spec- 
ified the profiles v(r), p(r) and AmiM for each stellar sub- 
component in each dark matter halo, we use Equation 1 1 of 
Merritt ( 1985) to calculate the corresponding phase-space dis- 
tribution functions. We check this calculation by performing 
N-body simulations in which stars orbit within the adopted 
potential and have initial positions/velocities drawn from the 
calculated distribution function. These simulations show no 
significant departures from the initial dynamical configuration 
after 100 crossing times, indicating that the calculated distri- 
bution functions indeed correspond to equilibrium dynamical 
models. 

Table [3] lists the grid of input parameters that specifies 60 
unique dynamical models that we use to represent individ- 
ual dSph stellar subcomponents. The top panel of Figure [3] 
displays the differential spatial distributions of stars and il- 
lustrates the effect of varying the slopes of the outer stellar 
density profiles over the range /3* =4,5,6. We note that the 
Plummer surface brightness pro files that we assume in our 
likelihood function (Section l3.2b correspond to (deprojected) 
stellar density profiles with (a*, 7*) = (2,5,0). Test cases 
with /3* = 4 and /3* = 6 correspond to stellar density profiles 
that decline less and more steeply, respectively, thereby vio- 
lating the simplistic assumption of Plummer surface bright- 
ness profiles that is inherent in our method. The bottom five 
panels in Figure [3] display the projected velocity dispersion 
profiles that we obtain from large {N = 10 5 ) random sam- 
ples of the corresponding distribution functions. These test 
models include velocity dispersion profiles that are flat, rising 
and/or falling, thereby violating the simplistic assumption of 
constant velocity dispersion that is inherent in our method. 

4.2. The Meaning of T 

For the cored and cusped dark matter halos that we con- 
sider in our tests, the top panel of Figure |4]displays enclosed- 
mass profiles M(r) and the bottom panel displays logarith- 
mic slopes dlogM /dlogr. This figure illustrates the sim- 



TABLE 3 

Tests on synthetic data: grid of input parameters for 
dynamical test models 



Profile 


Parameter 


values considered 


Stellar Subcomponent (Eq.|15t 








r*/r D M 


0.10,0.25,0.50,1.0,1.5 




a* 


2 




/3* 


4,5,6 




7* 


0.1 


Dark Matter Halo (Eq. [16} 


r a /r, 


1, 00 








Po/(M©pc- 3 ) 


0.064 




iDM/kpc 


1 




«DM 


1 






3 




TDM 


0, 1 




FlG. 3. — Top: Differential distribution of stars as a function of (projected) radius 
for the stellar subcomponents in our tests. Bottom five panels: Projected velocity dis- 
persion profiles for the physical dynamical models used to construct the synthetic data 
sets on which we test our method. For clarity we plot profiles corresponding only to 
models with outer stellar density profiles specified by /3* =5 (velocity dispersion pro- 
files for models with /3* = 4 and /3* = 6 behave similarly). Notice that for a given halo 
potential, projected velocity dispersion profiles corresponding to isotropic (r a = oo) and 
anisotropic (r a = r*) velocity distributions cross at r ~ r*, the radius where the number 
of stars reaches a maximum. This phenomenon helps to explain the insensitivity of our 
mass estimates to velocity anisotropy (Section [4.4. \\ . 
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FIG. 4. — Enclosed mass profiles (top) and slopes of logarithmic mass profiles (bot- 
tom), for the cored and cusped dark matter halos considered in tests of our method. 
Discrete points identify the luminous scale radii of various dSph-like stellar subcompo- 
nents that we embed in these halos for our tests. 

pie relationships between the inner slope of the logarith- 
mic density profile (7dm), the inner slope of the logarith- 
mic mass profile (lim r ^o[^logM/c/logr]), and the slope (T = 
AlogM/Alogr) that we actually measure. For the spher- 
ically symmetric dark matter density profiles specified by 
EquationfTBI the inner mass profile is given by lim, ^oM(r) oc 
r 3-7DM. Thus the inner slope of the logarithmic mass profile is 
given by 

lim[cflogM/<ilogr] = 3-7dm- (17) 

At their centers, cored (7dm = 0) dark matter halos have 
linv^ot^logM/aflogr] = 3 and cuspy NFW (7dm = 1) halos 
have lim,.^o[(f logM/iilogr] = 2. 

However, since our method constrains enclosed masses at 
the halflight radii of two stellar subcomponents, the slope 
T = A log M/ A log r that we measure corresponds to the 
slope at some finite radius r > and will not necessarily 
represent the central value of d log M/d log r. The bottom 
panel of Figure |4] illustrates that as r increases, the slope 
c/logM/c/logr decreases monotonically. Therefore an unbi- 
ased estimate of the slope T = AlogM/Alogr evaluated at 
radii r^z > Om > will necessarily be smaller than the central 
value of (ilogM/t/logr. In other words, an unbiased estimate 
of r constrains the central value of the slope to be 

\im[d\ogM/d logr] > T, (18) 

r— s-0 

and therefore, via Equation \T7\ constrains the inner slope of 
the logarithmic density profile to be 

7DM <3-r. (19) 

4.3. Synthetic Data 

Each of the 60 unique dynamical models included in the 
grid outlined in Table [3] describes a single equilibrium stel- 
lar component embedded in the potential generated by a dark 
matter halo. Therefore each model can correspond to either 
stellar subcomponent in a system that has two such subcom- 
ponents. We consider all possible combinations of scale radii 
(r*,i/r D M, r*,2/>"DM), outer slopes (/3*.i, /3*,2), and anisotropy 



radii (r a i, r a 2) for which both stellar subcomponents are 
embedded in the same dark matter halo. These combina- 
tions yield a total of 1080 unique, two-subcomponent struc- 
tural/dynamical models with which we test our method. In 
all cases where the two stellar subcomponents have differ- 
ent scale radii, we assign (see below) reduced magnesium 
indices such that the subcomponent with smaller r„ is the 
metal-rich subcomponent, consistent wi th the phenomenol- 
ogy of well-stu died Local Group dSphs (iTolstov et al.ll2004t 
Battag liaet all 120061: iBattaglia et al.1 1201 11) . When the two 
stellar subcomponents have the same scale radii (note that in 
these cases the slope V is undefined), we arbitrarily choose 
one to be metal-rich and the other to be metal-poor. 

For each of the 1080 unique structural/dynamical test mod- 
els we perform ten realizations. In setting up a given real- 
ization we draw stellar population parameters randomly from 
uniform distributions within the following limits: 

• sample sizes 3 < log I0 |7Vi + N2 + Nmw] < 4 (similar to 
the available MMFS samples) 

• member fractions 0.4 < (N\ +N2)/(N\ +N2+N M w) < 
0.9 

• subcomponent fractions 0.1 < Ni/(N\ +N2) < 0.9 

• mean systemic velocities (heliocentric rest frame) < 
(VO/Ckms- 1 ) <250 

• mean spectral index 0.3 < (W')i/A< 0.5 for the 
'metal-rich' subcomponent 

• mean spectral index separation < ((W')\ — 
(W') 2 )/k< 0.25 

• proper motions -100 < /i Q /(mas/cent) < +100 and 
-100 < ^/(mas/cent) < +100. 

We place half (randomly selected) of the synthetic 'dSphs' at 
the (3D) position of Fornax and the other half at the location 
of Sculptor (TableCD- 

With the above stellar population parameters specified for 
a given realization, we then use an accept/reject algorithm 
to draw the appropriate numbers of positions and velocities 
from discrete random samplings of the appropriate 6D distri- 
bution function. We then project the positions and velocities 
along the line of sight in order to mimic observables R and 
V. Next we assign reduced Mg indices, W', to each star ac- 
cording to whether it is drawn from the 'metal-rich' or 'metal- 
poor' subcomponent. We assign W' values to the metal-rich 
and metal-poor member stars by drawing values from Gaus- 
sian distributions with variances er^, j = o^, 2 = 0.02 A 2 and 
means drawn randomly from the ranges specified above. To 
the line-of-sight velocities of all member stars we apply red- 
shifts appropriate to the synthetic dSph's systemic 
3D space motion and line of sight (Equation [TO}. Finally, 
we scatter all velocities and W' values according to actual 
measurement errors drawn randomly from the MMFS data set 
(median errors are ty ~ 2 km s" 1 and ey?i ~ 0.01 A). 

To stars drawn from a 'foreground' contamination compo- 
nent we assign positions drawn randomly from a uniform spa- 
tial distribution (within the projected position of the outermost 
member star) and assign velocities drawn randomly from the 
Besangon model of Milky Way stars (filtered by our photo- 
metric criteria for selecting dSph red giants) along the line of 
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sight to the either (chosen randomly in each realization) the 
Fornax or the Sculptor dSph. To foreground stars we assign 
W' values and associated errors drawn directly from measure- 
ments of probable (P mem < 0. 1) foreground stars in the MMFS 
data set. 

4.4. Systematic Errors 

We apply our method (including the initial estimation of 
membersh ip probabilities using the single- com ponent EM al- 
gorithm of Wal ker et al.l l2009b: see section [272b to each of the 
10800 synthetic data sets (ten realizations for each of the 1080 
unique structural/dynamical models). For our tests we are less 
interested in distributions of parameter estimates than we are 
interested in distributions of errors E(S) = S-Si npnt . Since it is 
impractical to examine the error distribution obtained in each 
individual realization, in what follows we consider error dis- 
tributions obtained after 'stacking' (accomplished by drawing 
a fixed number of points from the final chains obtained for 
each of the individual realizations that are to be combined) er- 
ror distributions corresponding to input models that have var- 
ious parameters of interest (e.g., 7dm, fa/ f*, '"♦/'"dm) in com- 
mon. By combining error distributions in this way, the result- 
ing 'composite' error distributions effectively are marginal- 
ized over other input parameters (e.g., sample sizes, mem- 
ber fractions, spectral-index parameters, etc.) that could vary 
from realization to realization. 

After stacking results according to whether the input dark 
matter halo is cored or cusped, Figure [5] displays error distri- 
butions for ten of the twelve free parameters involved in our 
likelihood function (Equations [6] and Q31 since the test mod- 
els do not have constant velocity dispersion, errors associated 
with free parameters u\ and a\ are poorly defined and there- 
fore not plotted). For both cored and cusped input halos, we 
recover unbiased estimates of most of the free parameters in 
our likelihood function, as indicated by the symmetry about 
zero of the composite error distributions shown in Figure [5] 

Notable exceptions appear in panels that display error dis- 
tributions Efax/rhp) and E(log l0 [rf,3./pc]) for the parame- 
ters that specify halflight radii of the two stellar subcompo- 
nents. In general we over-/under-estimate halflight radii of 
subcomponents for which density profiles decline less/more 
steeply than the Plummer profiles (/3* = 5) adopted in our like- 
lihood function. For the outer subcomponent, the three peaks 
at ^(logjofr/j^/pc]) ~ -0.2, and +0.3 correspond to tests 
in which the input model has /3* .2 = 6, 5 and 4, respectively. 
Notice that the relative offsets among these peaks directly cor- 
respond to offsets in the radii at which the differential spatial 
distributions dN/dR are maximized for the three stellar den- 
sity profiles (Figure[3]i. The distribution of errors E(r/ U i / r/, 2) 
is centered at zero because we tend to recover accurate esti- 
mates of r/j i/r/, 2 for cases in which = /3».2, even when 
neither component declines like the assumed Plummer pro- 
file. When the stellar density profile of the inner subcompo- 
nent declines less (more) steeply than that of the outer sub- 
component, our method tends to over- (under-) estimate the 
ratio ri h \/ri h 2, thereby broadening the distribution of errors 
E( r h, iA/1,2) compared to what would be found if we consid- 
ered only cases with /3* 1 = /3*,2- Reassuringly, the central 
peaks at£'(log 10 [r/,.2/pc]) ~ and£'(r;,,i/r/ ! .2 ~ indicate that 
we recover unbiased estimates of halflight radii when the in- 
put model has the assumed Plummer values of /3* 1 = /3».2 = 5. 
It is also reassuring that even when our estimates of the sub- 
component halflight radii have significant errors, as they tend 



to when one or both of the stellar subcomponents in the input 
model have /3* =/5, we recover reasonably unbiased estimates 
of the other free parameters. 

4.4.1. Masses 

As mentioned above, errors associated with free parameters 
a\ and o\ are poorly defined because the physical dynamical 
models that we invoke in order to test our method generally 
do not have constant velocity dispersion. It is important to 
realize that this aspect of the test models violates not only the 
assumption of constant velocity dispersion (Equation |9]l that 
enters our likelihood function, but also the assumption of con- 
stant velocity dispersion that lea ds to the mass estim ator given 
by Equation |2] in the first place dWalker et al.ll2009ah . Before 
we examine the errors associated with the masses returned by 
our method, let us first examine the bias that is associated di- 
rectly with the mass estimator on which our method is based. 

For each of the 60 unique dynamical models (Table |3]i 
that we use to represent dSph stellar subcomponents, Figure 
[6] plots the difference between the true value of M(r/,) and 
the value we obtain by applying our simple mass estimator 
(Equation^ to an arbitrarily large — and therefore effectively 
noiseless — sample drawn randomly from the phase-space dis- 
tribution function. This figure establishes that there is bias 
associated directly with our simple mass estimator, and high- 
lights its important characteristics. 

First, the mass errors that are apparent in Figure [6] show 
no significant dependence on the velocity anisotropy of the 
input model. For a given dark matter halo and embedded- 
ness of the stellar component r*/roM, Equation [2] provides 
mass estimates having similar errors regardless of whether 
the tracer velocity distribution is isotropic (r a = 00) or has a 
radially variable anisotropy profile (r a = r»). This lack of a 
dependence on anisotropy can be understood phenomenolog- 
ically upon re-examination of Figure [3] There we find that 
for fixed 7dm and r*/ruM, projected velocity dispersion pro- 
files corresponding to isotropic and anisotropic distribution 
functions intersect at approximately the luminous scale ra- 
dius, R ~ r*. This is also the radius where the differential 
spatial distribution dN/dR reaches a maximum (top panel of 
Figure [3}. Therefore, samples drawn from both isotropic and 
the anisotropic distribution functions tend to be dominated by 
stars at radii where the difference in velocity dispersion be- 
tween the two cases is negligible. As a result, we measure 
similar velocity dispersions and therefore obtain similar mass 
estimates for isotropic and anisotropic cases. 

Second, Figure [6] shows that while the magnitude of the 
error depends on whether the dark matter halo is cored or 
cusped, for both types of halo the error varies with the em- 
beddedness of the stellar component (as quantified by the ra- 
tio r*/roM) in the same way. Specifically, values of M(r/,) ob- 
tained from Equation [2] are more strongly over-estimated for 
stellar components that are more deeply embedded (smaller 
r»/roM) in their respective dark matter halos. These tests in- 
dicate that this bias is inherent in mass estimators of the form 
Af(~ fh) oc rhO 2 (where a 2 is the global velocity dispersion). 

Having identified biases inherent in our mass estimator it- 
self, we now examine how those systematic errors propagate 
through our MCMC analysis. In order to see how our er- 
rors depend not only on the central slope of the dark mat- 
ter density profile, but also on velocity anisotropy and em- 
beddedness of the stellar subcomponent, we now stack er- 
ror distributions corresponding to input models that have the 
same values of 7dm, fa/ r * an d r t /r DM . We do not break 
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FIG. 5. — Recovery of the free parameters in our likelihood function (Equations[6]and [T4l . from tests with synthetic data. Panels display composite error distributions (red/blue for 
the distributions obtained by stacking error distributions obtained in individual realizations corresponding to cored/cusped input halos) evaluated by subtracting the known input value of 
each parameter from the MCMC-sampled values. Peaks at £(log 10 [r/,.2/pc]) ~ —0.2, 0.0, +0.3 correspond to input models with outer stellar density profiles specified by /3*,2 = 6, 5,4, 
respectively. Since the test models generally do not have constant velocity dispersion, errors associated with velocity dispersion estimates are poorly defined and therefore not shown 
(but see FigureQJfor errors associated with derived masses). 
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FIG. 6. — Errors associated with the simple mass estimator given by Equation[2] for 
the dynamical models that we use to construct synthetic data sets. Models are plotted 
according to whether the input halo is cored (7dm = 0) or cusped (7dm = 1), whether 
the input velocity distribution is isotropic (r a = oo) or has Osipkov-Merritt anisotropy 
(r a = r*), and whether the input stellar density profile has outer slope specified by 
/3* = 4,5,6. Regardless of the value of 7dm, ''a or /3* for the input dynamical model, 
errors £(log lo [M(r/,)/M ]) = log lo [M(r ( ,)/M0] es , imffled -log m [M(n,)/M e ] intmt corre- 
late primarily with the degree to which the stellar subcomponent is embedded within the 
dark matter halo. 



results down further according to the input values of /?* be- 
cause we find no significant dependence of the mass errors 
on this parameter (see, e.g., figure |6j. Figure [7] displays the 
stacked distributions of mass errors E(\og w [M(jh) / Mq\) = 
log 10 [M(r/,)/M Q ]-log 10 [M in p U t(r/ ! )/Mo]. Here we define the 
mass error as the difference between the estimated mass at the 
estimated half-light radius and the true mass at the estimated 
half-light radius, since for real dSphs one does not know the 
subcomponent halflight radii a priori. Figure|7]shows that we 
generally recover the systematic behavior of errors that we 
expect to result from the bias inherent in the mass estimator 
as previously exhibited in Figure [6] (indicated by downward- 
pointing arrows in Figure |7). More specifically, we recover 
the expected anti-correlation between £'(log 10 [M(r/,)]) and 
r*/ '"dm- 

We note that for our cored input models with r»/j"DM = 



0.1 and r*/roM = 0.25, the distributions of mass errors 
£'(log 10 [M(r/ ! )]) peak at values that are slightly smaller than 
the errors expected from the bias inherent in the mass esti- 
mator (top two rows of panels in Figure |7). This peculiarity 
follows from the fact that for these particular models, the in- 
put 'global' (i.e., calculated from a large sample of velocities 
drawn randomly from the subcomponent's distribution func- 
tion) velocity dispersions are ^2-4 km s" 1 (see Figure [3J, 
similar to the MMFS velocity errors adopted for our tests. 
We resolve these small dispersions only marginally and our 
resulting estimates are biased toward values that are smaller 
than the input global dispersions. The corresponding masses 
are therefore not overestimated as strongly as would be ex- 
pected from the bias that is inherent in the mass estimator. 
This complication does not affect the velocity dispersions (or 
masses) we measure for subcomponents in real dSphs (Sec- 
tion[5}, which are all > 6 km s" 1 and thus well resolved with 
MMFS data. 

4.4.2. Slope of the Mass Profile 

For the present study we are concerned primarily with the 
error associated with our measurement of the slope T. The 
most important result from our tests with synthetic data is 
the finding that regardless of whether the dark matter halo is 
cored or cusped, isotropic or anisotropic, our method tends 
to return more strongly over-estimated masses when the stel- 
lar subcomponent is more deeply embedded within the dark 
matter halo (Figures [6] and |7). Consequently, for a dSph 
with two stellar subcomponents having different scale radii, 
our method tends to overestimate the mass of the inner sub- 
component more strongly than it does the mass of the outer 
subcomponent. Thus we can expect that our estimate of the 
slope T = (log[M(r Ai2 )/Af(r Ai i)])/(log[n,,2/r A ,i]) will tend to 
be under-estimated. 

This is exactly what we find in our tests. Figure [8] 
displays^ composite distributions of slope errors E(T) = 
r-(log[Mi n p Ut (r/,,2)/Mi n p Ut (r/,. 1 )])/(log[r hi2 /r, li i]), stacked ac- 

8 In Figure [8] we exclude PDFs corresponding to structural/dynamical 
models in which r„ . i = r * %, since for these cases V is undefined. Accordingly, 
the PDFs for V obtained in these tests indicate no meaningful constraint. 
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FIG. 7. — Distributions of mass errors E(iog m [M(r/,)/MQ]) = log l0 [M(r/,)/MQ] - 
log !o [Mi npul (r/,)/M0 ] obtained for individual stellar subcomponents in tests of our 
method using synthetic data (Section l4.4.U . Here we have stacked error distributions 
from individual realizations corresponding to input models that have the same slope 
for the inner dark matter density profile (7dm)> anisotropy radius (r n ) and level of em- 
beddedness (r* //'dm) of the stellar subcomponent within the dark matter halo; we do 
not separate further according to the outer slope of the stellar density profiles be- 
cause we find no significant dependence of the mass errors on this quantity. Downward- 
pointing airows identify the mean error that is associated purely with the mass estimator 
itself (see Figure l6l. 



cording to the input value of 7dm- As expected, the anti- 
correlation between £'(log 10 [M(r/ 1 )]) and r^/r^M causes the 
bulk of these distributions to have values E(T) < 0, indicat- 
ing that indeed we tend to underestimate V for both cored and 
cusped input halos. 

We note that error distributions obtained in individual re- 
alizations of all tested models have finite widths and tails 
that include some positive errors E(T) > 0. However, of 
the models that we consider, only those with input values of 
7dm = and r*i/i"mA < 0.25 have error distributions that are 
sometimes centered — as quantified by either the maximum- 
or median-likelihood value — at E{T) > 0. These particular 
cases are responsible for generating the relatively thick tail 
toward positive values seen in Figure[8]for the distribution of 
E{T) corresponding to cored input models. As with the mass 
errors discussed in the previous section, these peculiarities re- 
sult from sampling rather than systematic error, since for these 
particular models both subcomponent velocity dispersions are 
small (~ 2-4 km s -1 , see Figure[3]l and only marginally re- 
solved with MMFS-quality data. Again, this resolution is- 
sue does not affect our analysis of real dSphs (Section |5), for 
which we measure subcomponent velocity dispersions of > 6 
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FIG. 8. — Top-Right: Distributions of slope errors E{Y) obtained in tests of our 
method using synthetic data (Section |4.4. 2k Here we have stacked error distributions 
from individual realizations corresponding to input models that have the same slopes for 
the inner dark matter density profile (7dm). We do not include results for test cases in 
which the input stellar components have the same halflight radii, since the slope V is 
then undefined. In general we identify a bias such that our estimates of the logarithmic 
slope of the mass profile generally are smaller than known input values (Section [4.4.2l , 
as expected from the dependence on embeddedness r* /rDM of the bias associated with 
our mass estimator (Figure |6|. 

km s -1 , at least three times larger than the median MMFS ve- 
locity measurement error. 

4.5. Summary of Test Results 

We summarize the results of our tests with the following 
conclusions. For the ranges of dynamical models and stellar 
population parameters considered here: 

• our method provides reasonably accurate and unbi- 
ased estimates of stellar population parameters f r 



/sub, (W')u(W') 
(Figure 0. 



mem > 

2 



-(W)2,log 10 [c r w',l] andl °gl0[c r W',2] 



• Our method tends to over-/under-estimate halflight radii 
of stellar subcomponents that have density profiles de- 
clining less/more steeply than the Plummer profiles as- 
sumed in our likelihood function (Figure |5j. 

• The assumption of constant velocity dispersion that is 
inherent in our mass estimator (Equation |2]i results in 
bias in our estimates of the masses M(r/,) of both stel- 
lar subcomponents (Figure[6]). This bias shows no sig- 
nificant correlation with the velocity anisotropy of the 
stellar subcomponents (Figures [6] and [7]). Rather, this 
bias correlates primarily with the degree to which the 
stellar subcomponent is embedded in the dark matter 
halo, such that (so long as both velocity dispersions are 
resolved) we tend to overestimate the mass enclosed 
within the halflight radius of the inner subcomponent 
more strongly than we do the mass enclosed within the 
halflight radius of the outer subcomponent (Figures [6] 
and|7|. 



This anti-correlation between mass error 
(£'(log 10 [M(r/ 1 )/MQ])) and degree of embedded- 
ness (V*/Vdm) of the stellar subcomponent results in 
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a tendency to Mnt/er-estimate the slope of the mass 
profile for both cored and cusped input halos (Figure 

Furthermore, we have already seen (Section l4~2l that the value 
of the slope T = AlogM/Alogr evaluated at two points 
Oi,2 > r h,\ > is smaller than the central value of the slope 
of the logarithmic mass profile (Figure HJ. The fact that we 
tend to underestimate V implies that Inequalities [18] and [19] 
continue to hold if the true value of T is replaced by our esti- 
mate. 

We note that the tests described in this section cover a fi- 
nite range of dynamical models. In principle there is no limit 
to the number and variety of models that might be applied to 
dSphs. In addition to the models described in this section, we 
have performed a similar set of tests using dynamical models 
that allow for either constant radial (/? an ; = +0.25) or constant 
tangential (/3 an j = -0.45) velocity anisotropy. The correspond- 
ing distribution functions were sampled b y Mark Wilkinson 
(privat e communication), as described by Charbonni er et al.l 
(1201 11) . In response to a request from the referee, we have also 
tested a model that reproduces the velocity dispersion profiles 
(particularly the st eeply falling p r ofile of the inner compo- 
nent) measured by Battaglia et al. (2008a) for metal-rich and 
metal-poor subcomponents in the Sculptor dSph. We find that 
we can reproduce these profiles with a model in which in- 
ner (outer) subcomponents follow King] (119621) surface den- 
sity profiles with core radii r c = 195 pc (r c = 250 pc), tidal radii 
r t = 6r c (r, = 15r c ), Osipkov-Merritt anisotropy radii r a = 195 
pc (r a = 500 pc), and are embedded in a cored (7dm = 0) dark 
matter halo with central density po = 0.5M Q /pc 3 and scale 
radius /dm = 360 pc. In our tests of these additional dynam- 
ical models we continue to find the same tendency to over- 
estimate M(r/,) according to the degree to which the stellar 
component is embedded within the dark matter halo, and thus 
the same tendency to underestimate V. 

Finally, we note that while the test models that we have 
considered include complexities (rising/falling velocity dis- 
persion profiles and non-Plummer stellar density profiles) that 
violate the assumptions of our method and give rise to the 
systematic errors that we have identified, there is no doubt 
that real dSphs have additional complexities that we have not 
simulated in our tests. Examples may include non-sphericity 
of stellar and dark matter distributions, non-equilibrium kine- 
matics, internal rotational components, binary stars, the pres- 
ence of three or more distinct stellar subcomponents, etc. In 
Section[6]we discuss the potential for sensitivity of our results 
to these possibilities. 

5. RESULTS FOR CARINA, FORNAX AND SCULPTOR 

We now apply our method to the published MMFS data for 
the Carina, Fornax and Sculptor dSphs. The top three rows 
of panels in Figure [9] display posterior PDFs for each model 
parameter. Table |4] lists for each parameter the median value 
from the posterior PDF, with error bars indicating intervals 
that enclose the central 68% (and 95%) of values. 

5.1. Carina 

Our method returns no compelling evidence for chemo- 
dynamically distinct stellar subcomponents in Carina, for 
which we obtain / su b = 0.04^qq[ (the 95% error is consistent 
with / su b = 0). This non-detection may be surprising given 
Carina's episodic star fo rmation history as eviden ced by its 
two horizontal branches ( Smecker-Hane et al.l 19941) and three 



main-sequence turn-offs dHurlev-Keller et aTl 11998b . How- 
ever, a previous analysis bv lKoch et al.l (I2006T) of an indepen- 
dent VLT/FLAMES data set uncovers only a weak metallic - 
ity gradient i n Carina, a result qualitatively repr oduced using 
MMFS data dWalker. Mateo & Olszewskill2009l) . Thus while 
Carina may in fact have independent stellar populations, the 
apparently weak coupling of chemical and dynamical proper- 
ties prevents a clear separation using our method. As we do 
not clearly separate two stellar subcomponents, we obtain no 
meaningful estimate of T for Carina. 

For the single stellar component that includes a fraction 
/mem = 0.41*oo[ of stars in the Carina field, our method returns 
estimates of the halflight radius (r/, = 260*.}q pc) and veloc- 
ity dispersion (ay = 6.4*[} j km s" 1 ) that agree with previously 
published values based on analyses that assume a single stellar 
component (e.g.. iMateo et al.l Il993t llrwin & Hatzidimitrioul 
119951: iMunoz et al.ll2006l) . The proper motion returned by our 
method (/j, n = +40*j}7 mas /century, fig = +17*15 mas/century; 
see also lWalker et al]l2008l) is weakly c onstrained but agre es 
with the HST astrometric measurement (Pia tek et al.l l2003). 

5.2. Fornax and Sculptor 

For the remainder of this section we shall consider only For- 
nax and Sculptor, for both of which our method clearly dis- 
tinguishes two stellar subcomponents. We find that a frac- 
tion / su b = 0.60±g ■['[ of Fornax members belong to a rela- 
tively metal-rich subcomponent ((W')\ = 0.47*y{j[ A) with 
r/,j = 559*32 P c an d velocity dispersion oy.i = 10. O*]}^ km s _1 , 
while the metal-poor ((W'} 2 =031^1 A) subcomponent has 
r h,2 = 888*| j pc and velocity dispersion <r V 2 = 14.5*2 S km s -1 . 
For Sculptor we find that a fraction / su b = 0.53*Q g of mem- 
bers belong to a metal-rich ((W')\ = 0.36*q qJ A) subcompo- 
nent with r/,.i = 167*j pc and velocity dispersion oy.i = 6. 5*^5 
km s" 1 , while the remaining metal-poor ((W'}2 = 0.28*q{]} 
A) subcomponent has r/, 2 = 302*^ pc and velocity disper- 
sion oy ? = 1 1 .6*°^ km s " 1 . Thus we re c over t he results of 
iTolstoy et all (120041) and IBattaglia etakl (120061) . who show 
that the metal-rich inner subcomponents of Sculptor and For- 
nax have smaller velocity dispersions than the metal-poor 
outer subcomponents (furt hermore, the sc a le radii of 150 ± 20 
pc and 350 ± 10 pc that IBattaglia etail J2008a) derive for 
inner and outer scale radii, respectively, are broadly consis- 
tent with our measurements). The constraints on Fornax's 
proper motion stand in excellent agreement w ith indepen- 
dent a strometric measurements made with HST (Piate k et alj 
120071) . Our measurement of Sculptor's pr oper motion dis- 
agrees wit h astrometric mea surements by Schweitz er et alj 
( 1 19951) and[Piatek et al. (2006), which also disagree with each 
other. If either of the astrometric measurements is correct, 
then the velocity gradient that is pre sent in Sculptor may in- 
dicate a small rotational component (Battaglia et al. 2008a). 

For both galaxies the estimates of halflight radii and veloc- 
ity dispersions simultaneously provide estimates of masses 
\og w [Mi(r hA )/pc], log 10 [M2(r/,.2)/pc] (Equation O, and the 
slope r = A logM/ A log r (Equation[5]). Left and center pan- 
els in Figure [10] display our MCMC-sampled values in the 
plane of halflight radius and enclosed mass, color-coded ac- 
cording to the likelihood obtained from Equation [14] (nor- 
malized by the maximum-likelihood value). The right-most 
panel of Figure[l0]displays the corresponding posterior PDFs 
for r for both Fornax and Sculptor. Formally, we obtain 
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T = 2.6V$$ for Fornax and V = 2.95^ for Sculptor. 

5.3. Significance 

Recall from Section 14.21 that because we estimate masses 
at two finite points r/,,2 > rh.\ > 0, the resulting slope T = 
AlogM/Alogr corresponds to a region where the instanta- 
neous slope d\ogMld\ogr is smaller than the central value 
of 3 - 7dm (Figure [4j- Furthermore, the amount by which 
it is smaller depends on the scale radius of the dark matter 
halo, a quantity that we do not attempt to constrain here. As 
a result, a particular value of T is strictly inconsistent only 
with central logarithmic density slopes 7dm that are larger 
than 3 — r (Inequality [T9l. This means that measured values 
in the range 2 < T < 3 are consistent with cores (7dm = 0; 
c/logM/c/logr < 3 at all nonzero radii), but inconsistent with 
NFW and steeper cusps (7dm > l;c/logM/c/logr < 2) at all 
nonzero radii). 

On these grounds we can use the posterior PDFs for T, de- 
noted P(T), to calculate the significance, s(7dm), with which 
our measurements exclude dark matter density profiles with 
central slopes equal to or steeper than a particular value of 

7dm: 

/oo 
P(T)dT 

*(7dm)=1— ^ ■ (20) 

/ P(T)dT 

J —oo 

Our measurements exclude NFW and steeper cusps at 
s(7dm > 1) > 95.9% (Fornax) and s( 7dm > 1) > 99.8% 
(Sculptor) significance levels. For several reasons we regard 
these formal exclusion levels as conservative. First, for the 
equilibrium dynamical systems that we consider as test cases 
in Section @] our method tends to overestimate the mass of 
the inner subcomponent more strongly than it ove restimates 
the mass of the outer subcompone nt (Sec tion l4.4.T1 i. resulting 
in an underestimate of T (Section [4.4.2b . The only counter- 
examples to this trend that we find in our te sts res u lt from 
poorly resolved velocity dispersions (Sections 14.4. 1 1 - [4.4. 2\ . 
This resolution problem does not affect our results for For- 
nax and Sculptor, where even the smallest velocity dispersion 
that we measure (ay i = 6.5±Q5 for Sculptor's inner subcom- 
ponent) is several times larger than the median velocity error 
in the MMFS data set. 

Second, by setting the lower limit of the integra- 
tion in Equation [20] at the value of the central slope 
lim,._j.()[a! log M/d \ogr] = 2, we extend maximum generosity 
to models with 7dm > 1, which have instan tane ous slopes 
lilogM/iilogr < 2 at all nonzero radii (Section l4~2l and Figure 
|4j. At the radii (> 300pc) where we evaluate V for Fornax and 
Sculptor, the highest-resolution Aquarius simulations predict 
that dSph-like CDM halos have dlogp/d log r ~ 1.3, orequiv- 
alentl y, lilogM/iilogr ~ 1.7 (see Figure 23 of ISpringel et al.l 
2008). Our measurements rule out these slopes with signifi- 
cance > 99.54% (Fornax) and > 99.97% (Sculptor). 

Third, we have assumed that the stellar subcomponents 
contribute negligibly to the gravitational potential. This as- 
sumption generally holds for dSphs, but least so for For- 
nax, where the dynamical m ass-to-light ratio is M/Ly ~ 10 
in solar units (lMateol[l998l) . If we attempt to remove the 
stellar contribution to the enclosed mass at each radius us- 
ing the best-fit Plummer profiles to both stellar subcompo- 
nents, we find that for any plausible stellar mass-to-light ratio 
0.5 < M/L V /[M/L V ] Q < 5, our estimates of V increase by 



a few percent (because the stars contribute a larger fraction 
of mass to the inner than to the outer point), again exacer- 
bating the discrepancy with halo models having 7dm > 1 . In 
summary, all systematic errors that we have identified behave 
such that the significance levels we report are conservative. 

Finally, we note that for dark matter density profiles of the 
form given by Equation [T6l values of dlogM /dlogr > 3 are 
unphysical, as they imply 7dm < (Inequality [T9|>. We note 
that our method does not rule out such unphysical values, 
which is unsurprising since we have not imposed any phys- 
icality constraints. However, it is reassuring that the bulk of 
our posterior PDFs correspond to physically plausible scenar- 
ios with T < 3. 

6. DISCUSSION 

Let us review the assumptions that enter into our measure- 
ment of T. In formulating our method we assume that a dSph 
consists of either one or two spherically symmetric, equilib- 
rium stellar subcomponents that independently trace the same 
spherical dark matter potential. In order to quantify proba- 
bility distributions for observed quantities, we further assume 
that both stellar subcomponents have Plummer surface bright- 
ness profiles, Gaussian Mg-index distributions, and Gaussian 
line-of-sight velocity distributions with constant dispersions 
that receive negligible contributions from 'non-thermal' phe- 
nomena such as rotational support and/or binary-orbital mo- 
tions. The tests described in Section[4]indicate that for a range 
of models that explicitly violate our assumptions about Plum- 
mer surface brightness and constant velocity dispersion pro- 
files, our method tends to underestimate T, implying that the 
stated NFW exclusion limits are conservative. Here we dis- 
cuss the potential for sensitivity to several assumptions inher- 
ent in our method that are not violated in the tests of Section 
H]but might be violated by real dSphs. 

6.1. Spherical symmetry 

Fornax and Sculpt or both have projected minor - to-major 
axis ratios of ~ 0.7 (llrwin & Hatzidimitrioul [l995h and are 
among the roundest of the Milky Way's dSph satellites. In 
order to investigate the degree to which the observed flatten- 
ing of Fornax and Sculptor might affect our measurements of 
r, we repeated our analysis using elliptical instead of circular 
radii, where a star's 'elliptical radius' is the semi-major axis of 
the ellipse (with center list ed in Table Q] position angle an d el- 
lipticity listed in Table 2 of llrwin &^ HatzidimitriouTl 995l) that 
passes through the position of the star. Use of elliptical instead 
of circular radii gives constraints of T = 2.72^{J 43 for Fornax 
(exclusion significance s(7dm > 1) > 96.1%) and T = 2 .40^ 
for Sculptor (exclusion significance s(7dm > 1) > 93.9%). 
Thus the NFW exclusion level for Fornax is relatively robust 
while the exclusion level for Sculptor shows mild sensitivity 
to whether or not we adjust for Sculptor's elliptical morphol- 
ogy- 

6.2. Dynamic equilibrium 

We measure V by twice applying a mass estimator (Equa- 
tion [2]l derived from the spherical Jeans equation (Equation 
Q}, which holds for spherical systems in dynamic equilib- 
rium. We must therefore consider the fact that nearby dSphs 
are vulnerable to disruptive tidal forces as they orbit within 
the potential of the Milky Way. Tides can temporarily inflate 
dSph velocity dispersions — and thus dynamical masses — by 
'heating' the stellar component during a close pericentric pas- 
sage and/or by stripping stars that, thence unbound, linger 
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FIG. 9. — Results for the Carina, Fornax and Sculptor dSphs. Panels display posterior PDFs for model parameters, obtained from applying the two stellar subcomponent models 
introduced in Section^] Table[2]lists median values and 68% (95%) confidence intervals derived from these PDFs. 
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FIG. 10. — Left, center: Constraints on halflight radii and masses enclosed therein, for two independent stellar subcomponents in the Fornax and Sculptor dSphs. Plotted points 
come directly from our final MCMC chains, and color indicates relative likelihood (normalized by the maximum-likelihood value). Overplotted are straight lines indicating the central 
(and therefore maximum) slopes of cored (linv-^o^logM/dlogr] = 3) and cusped {\\m r ^od\ogM j d\ogr\ = 2) dark matter halos. Right: Posterior PDFs for the slope Y obtained for 
Fornax and Sculptor. The vertical dotted line marks the maximum (i.e., central) value of an NFW profile (i.e., cusp with 7dm = 1, linv—^oldlogM/dlogr] = 2). These measurements 
rule out NFW and/or steeper cusps (7dm > 1) with significance 5 > 96% (Fornax) and s > 99% (Sculptor). 



sufficiently near the dSph to be observed and counted as 
bound members (e.g.. IPiatek & Prvorill995t lOhet al.lll995b 



Read et al. 2006; Kli mentowski et alJ 12007: Penarrubi a et al.1 
2008b, 2009). Both phenomena affect the outer more than 
the inner parts of a satellite — thus tidal heating is the only 
process we identify that may cause our method to return an 
over-estimate of T. 

However, measurements of their systemic dista nces and ve- 
locities imply that neither Fornax (D ~ 138 kpc, Mate d 19981) 
nor Sculptor (D ~ 79 kpc) experience strong tidal encoun- 
ters with the Mi lky Way. Fornax's line-of-sight velocity and 
proper motion ( Piatek et al] 120071 supported by this work) 
imply a pericenter distance of r p = 1 1 8+52 kpc (IPiatek et alJ 
2007, error bars give 95% confidence intervals), and Sculp- 
tor's imply r p ~ 65 kpc (with 95% confidence intervals al- 



lowing values as low as ~ 30 kpc) fo r either of the two astro- 
metric proper m otion measurements (Schweitzer et al.ll 1 995b 
Piatek et al ]l2006h. N-body simul ations by Penarrub ia et alJ 
(2009) and Penarrub ia et aLl dTolO) demonstrate that for satel- 
lite halos that follow the generic density profile given by 
Equation Q21 the instantaneous tidal radius at pericenter is 
r, w r ; ,[M dsph (< r,)/(3M MW (< >>)] 1/3 , where M dsph (r f ) is the 
dSph mass enclosed within the tidal radius and Mmw(< r p ) 
is the enclosed mass of the M i lky W ay within the peri- 
centric distance. Watki ns et alJ ([2010) have recently used 
a sample of tracers (halo stars, globular clusters and satel- 
lite galaxies) in the outer Galactic halo to estimate a mass 
of M MW (< 300kpc) = 0.9 ± 0.3 x 1O 12 M . We obtain con- 
servative lower limits for the pericentric tidal radii of For- 
nax and Sculptor by considering only the stellar mass of 
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TABLE 4 

Constraints from Magellan/MMFS Data: Carina, Fornax and Sculptor* 



Carina 



Fornax 



Sculptor 



Model Parameters 



Derived Quantities 



/mem 
/sub 

1i.l/>n,2 

logiolnu/pc] 

Wl/A 

({W) ] -{W) 2 )/A 

lo glo[°iu/A 2 ] 
togiot ^/^ 2 ] 

logiotoivo™ 12 ^ 2 )] 

logKjlo^AkmV 2 )] 
\x a j (mas / century) 
fj,g I (mas / century) 

log lo [M(r h ,i)/M ] 
log lo [A/(rh >2 )/M0] 
T = AlogM/Alogr 

3-r" 



0.41 



0.04 



0.83 



+0.0K+0.03) 
-0.01(-0.03) 

+0.01(+0.03) 
■0.01(-0.04) 

+0.12(40.17) 
-0.16(-0.36) 



9 47 +0.02(+0.03) 
Z ' q ' Z -0.02(-0.03) 

n 4 -+0.09(+O.20) 
"■^*-0.08(-0.13) 

n M +0.09(+0.20) 
u - -0.08(-0.13) 

_-, . c+I.10(+1.68) 
i l -1.25(-1.75) 



-2.48 



+0.10(40.18) 
-0.12(-0.29) 



4.46 



+0. 19(+0.38) 
O.28(-5.00) 



1.62 



+0.03(+0.07) 
-0.03(-0.06) 



40" 



+47(+92) 
-46(-92) 



17" 



+35(+70) 
-36(-69) 



6.97 



+0.04(+0.07) 
-0.04(-0.07) 



0.96 
0.60 



0.62 



2.95 



+0.01(+0.01) 
-O.Ol(-O.Ol) 

,+0.10(40.18) 
0.11 (-0.25) 

+0.05(4-0.09) 
-0.04(-O.09) 

+0.04(+0.08) 
O.03(-0.06) 



n / , 7 +0.01(+0.02) 
"•^'-0.01(-0.02) 

n , fi +0.02(40.05) 
-0.02(-0.05) 



-2.27" 



-2.07; 

2.00 

2.32 



+0.07(+0.14) 
-0.O9(-O.33) 

+0.12(+0.24) 
O.19(-0.63) 

,+0.05(+0.09) 
O.06(-0.16) 

+0.04(+0.08) 



0.04(-0.07) 



fi ,+14(+27) 
OJ -14(-28) 

_ 7 n+10(+20) 
-10(-21) 



7.67 



;.20 



2.61 



0.39 



+0.07(+0.12) 
-O.08(-0.20) 

+0.06(+0.13) 
-O.06(-0.12) 

+0.43(+1.07) 
-0.37(-O.68) 

,+0.37(+0.68) 
-0.43(-1.07) 




0.53 
0.55 



+0.0K+0.02) 
-O.01(-0.02) 

+0.07(40.14) 
-O.08(-0.15) 

+0.06(+0.12) 
-0.05(-0.10) 



9 .o+0.04(+0.09) 
Z - HO -0.03(-0.06) 

ft , fi +0.01(+0.02) 
u - JU -0.01(-0.02) 

n n „+0.01(+0.02) 
u ' uo -0.01(-0.02) 



-2.23 



+0.09(+0.18) 
-O.ll(-0.23) 



-3.30 



| +0.37(+0.56) 
'-0.90(-1.58) 



1.62 



40.06(40.11) 
-0.06(-0.13) 



2.13 



+0.05(+0.10) 
-0.04(-0.08) 



-33; 



+-26(4-51) 
-27(-54) 



_ 79 +34(+67) 
' -33(-67) 



6.77 



7.53 



2.95 



0.05 



+0.07(40.13) 
-O.07(-0.15) 

+O.08(+0.17) 
-O.07(-0.13) 

+0.51(+1.22) 
-O.39(-0.70) 

40.39(+0.70) 
-0.51(-1.22) 



Error bars enclose the central 68% (95%) of area under the marginalized ID posterior probability distribution function. 
From Inequality ! 19| this quantity represents an upper limit on the inner slope 7dm of the logarithmic density profile. 



each dSph (dynamical masses are > 10 times larger, Mateo 
[19981) and taking M M w(< 300kpc) as an upper limit for the 
mass enclosed within the pericentric radius of each satel- 
lite: r, > rp[L dsph T»/(3M M w(< SOOkpc))] 1 / 3 , where T» is 
the stellar mass-to-light ratio. Using the pericentric distances 
quoted above, the luminosities listed in Table[T]and assuming 
T* = IMq /Lq, we obtain r, > 2 kpc for Fornax and r, > 500 
pc for Sculptor. These lower limits are larger than the halflight 
radii that we measure for the outer subcomponents of both 
galaxies. Nearly all member stars lie inside these radii and 
thus — particularly if both galaxies additionally have dark mat- 
ter halos — we can expect that tides do not profoundly alter the 
structure and kinematics of Fornax and Sculptor during their 
pericentric passages. 

Even if Fornax and Sculptor have orbital pericenters much 
smaller than indicated by their measured proper motions, both 
systems are now sufficiently far from the Milky Way that 
they will have reached new equilibrium configurations and 
shed any stripped stars. Pefiarrubia et al] (12009b demonstrate 
that unbound tidal debris lingers near a tidally stripped satel- 
lite only for a time similar to the satellite's internal cross- 
ing time, which for the Milky Way's 'classical' dSphs is 
t c ~ (300pc)/(10kms _I ) w 30 Myr, or the amount of time re- 
quired to travel 10 kpc at constant speed 300 km s" 1 . This 
duration is much smaller than the amount of time required 
by either Fornax or Sculptor to travel to their present loca- 



tions from pericentric distances compatible with significant 
tidal disruption. 

6.3. Rotation 

Mass estimates for stellar subcomponents identified by our 
method are directly proportional to the corresponding esti- 
mates of stellar velocity dispersions. In principle, any con- 
tribution to these velocity dispersions by 'non-thermal' mo- 
tions such as rotational support or unresolved binary orbital 
motions (next section) might introduce a bias in our mass esti- 
mates beyond those that we have already identified in Section 

A stellar subcomponent that receives significant support 
against gravity from rotation about an axis not aligned with 
the line of sight will exhibit a smooth variation in mean veloc- 
ity as a function of position. For the simplest (solid body) ro- 
tation models, rotation introduces a gradient in mean line-of- 
sight velocity. All three of the dSphs studied here exhibit sta- 
tistically significant gradients in their velocities as measured 
in the heliocentric and Milk y Way rest frames (IWalker et al.l 
2008; Battaglia et al. 2008a). However, our method attributes 
any such gradient not to rotation (which we implicitly assume 
is insignificant), but wholly to the perspective effect induced 
by the d Sph' s systemic motion transverse to the line of sight 
(Section [3.3b . Since we account for this effect in our likeli- 
hood function, our method effectively removes the contribu- 
tion of any apparent velocity gradient from our estimates of 
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the subcomponent velocity dispersions. 

One might object that such gradients can arise due to a com- 
bination of perspective effects and rotation, and that by at- 
tributing any detected gradients entirely to perspective effects, 
we unduly ignore what might be real and dynamically signif- 
icant rotation. This concern is particularly relevant for Sculp- 
tor, where the proper motion that we estima te disagrees with 
both published astrome tric measurements dSchweitzer et alj 
ll995HPiatek et alj|2006l) . However, we find that the gradients 
in question are sufficiently small that our estimates of T are 
insensitive to our assumptions about the relative importance 
of rotation and perspective effects. If we assume that Fornax 
and Sculptor have zero proper motion in the heliocentric rest 
frame, such that any velocity gradient is allowed to contribute 
maximally to our velocity dispersion estimates, then we esti- 
mate T = 2.53^ for Fornax and T = 2.94^ for Sculptor; 
both values are statistically indistinguishable from our origi- 
nal constraints. 

6.4. Binary Stars 

It has long been recognized that dSph mass estimates de- 
rived from stellar velocity dispersions may be vulnerable to 
the in flation of those dispersi ons by stellar binary orbital m o- 
tions. lOlszewski et all d 19961) and lHargreaves et all d 19961) in- 
dependently demonstrate that unless dSphs have binary or- 
bital distributions that are pathologically different from those 
observed among field stars, binary motions do not contribute 
significantly to the velocity dispersions of ~ 10 km s" 1 mea- 
sured from red giants in 'clas sical' dSphs. On the other hand, 
iMcConnachie & Cotd (l2010h have recently shown that bina- 
ries may contribute significantly to (and in some cases dom- 
inate) the velocity dispersions of < 3 - 4 km s" 1 that have 
been measured from fainte r stars in some of the least luminous 
Milky Way satellites (e . g..lMartin et alj|2007t ISimon & Gehal 
2007t lAdenetai]l20 09: Kop osov et al.ll2011l). 

The Monte Carlo simulations of IMcConnachie & Cotd 
(2010) show that susceptibility to inflation by binary motions 
increases as the measured velocity dispersion decreases — 
i.e., for smaller dispersions it becomes more likely that bi- 
nary motions contribute significantly to the signal. Accord- 
ing to these simulations (and independent simulations by 
iMinor et"afl 120101) it is highly unlikely that binary motions 
contribute significantly to the velocity dispersions we esti- 
mate in this study, since the smallest dispersion we estimate is 
ay = 6.5^q 5 km s" 1 (for Sculptor's inner subcomponent). Fur- 
thermore, even if binary motions are significant in our sample, 
they would contribute more significantly to the smaller disper- 
sions that we measure for the inner subcomponents of Fornax 
and Sculptor than they do to the larger dispersions that we 
measure for outer subcomponents. Thus binaries would affect 
our measurement of the slope T in the same sense as t he bia s 
that is already inherent in our mass estimator (Section |4Aj}, 
causing us to overestimate the mass of the inner component 
to a larger degree than we overestimate the mass of the outer 
component. Therefore, to the extent that binary motions are 
relevant at all, they join the other sources of error that cause 
us to underestimate V systematically, again implying that our 
formal exclusion levels are conservative. 

6.5. Number of stellar subcomponents 

Our method allows for up to two stellar subcomponents. 
Deep photometric surveys suggest that Sculptor probably 
has exactly two stellar subcomponents that are sufficiently 



bright to observe in large numbers (Sculptor also has ex- 
tremely inetaFgoor stars similar to those found in the faintest 
dSph s dKirbv et al.ll2009t iFrebel et al.ll20Tot iTafelmever etaTI 
2010), which may correspond to a third, 'ultrafaint' subcom- 
ponent), as indicated by the different sp atial distributions of 
its blue- and red-horizontal- branch stars dHurley-K eller et al. 
1999; Majewski et al. 1999). For nax exhibits the same sort 
of bimodal horizontal branch morphology, but while Sculptor 
stopped forming stars as long as ~ 10 Gyr ago, Fornax has 
a you ng main sequence population with age o f order ^100 
Myr dStetson et al.lll998l iBattaglia et al.1 120061) . The young 
main sequence stars are more centrally concentrated than ei- 
ther of the subcomponents traced by evolved giants, so Fornax 
appears to have at lea st three distinct stellar subcomponents 
(BattagliaitiniSOOl). While the presence of three subcom- 
ponents offers the possibility of measuring the slope of For- 
nax's mass profile using three resolved points, it is unlikely 
that our spectroscopic sample of stars selected from Fornax's 
red giant branch includes a significant number from the young 
and relatively unevolved population. 

6.6. Future Work 

Finally, we note that the Local Group hosts several more 
dSphs for which similar measurements are feasible in the 
short term. Although its large angular size presents obser- 
vational challenges, there is already spectroscopic evidence 
that the Sextans dSph contains chemo-dyna mically distinct 
stellar subcomponents (Batta glia et al.l [201 1). Furthermore, 
star formation histories measured by comparing observed and 
synthetic color magnitude diagrams indicate that the Milky 
Way dSphs Draco, Leo I, Leo II and Ursa Minor all have in- 
termediate (1 < age /Gyr < 8) as well as old (age > 10 Gyr ) 
stellar populat ions (Dolp hin et aLll2005h ribistov et al.ll2009l) . 
lHarbeck et al.l (12001b identify spatially segregated blue- and 
red-horizontal branch p opulations in Andromed a satellites 
And I and And VI, and Mc Connachie et alJ (120071) find sim- 
ilar photometric evidence for distinct stellar subcomponents 
in And II. The ability to separate these objects statistically 
into chemo-dynamically distinct stellar subcomponents will 
depend on a combination of spectroscopic sample size and 
the individual properties of the dSphs, particularly regarding 
the inherent contrast between stellar subcomponents. A large 
and comprehensive survey of dSphs that are separable into 
stellar subcomponents will help to evaluate the generality of 
our results for Fornax and Sculptor. 

7. SUMMARY 

We have introduced a method for measuring the slopes of 
mass profiles for dSphs that have chemo-dynamically distinct 
stellar subcomponents. The method operates directly on spec- 
troscopic data and invokes neither a dark matter halo model 
nor any assumption about velocity anisotropy. Rather, it ex- 
ploits the basic principle that two resolved points in the same 
profile are sufficient to define a slope. Given measurements of 
stellar positions, velocities and spectral indices, our method 
can measure halflight radii and velocity dispersions simulta- 
neously for up to two stellar subcomponents in the same dSph. 
For any mass estimator of the form M(k^) oc ryjy, where k 
is some constant, these measurements immediately provide an 
estimate of the slope V = A log M / A log r as defined by two 
points. 

Applying our method to published Magellan/MMFS data, 
we distinguish two stellar subcomponents each in the Fornax 
and Sculptor dSphs, for which we measure T = 2. 61^ 37 and 
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r = 2.95^o 39 respectively. These values are consistent with 
cored dark matter halos of constant density over the central 
few hundred parsecs and rule out NFW-like (d log M/d log r < 
2) cusps with significance > 96% and > 99%, respectively. 
Tests with synthetic data drawn from physical distribution 
functions demonstrate that these exclusion levels are conser- 
vative. These results provide direct evidence against the no- 
tion that the current CDM paradigm successfully accounts for 
the phenomenology of dark matter at small scales. 
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